function h = CIC(w, R, N, D)
    //w -> Normalized frequncy point (1/Fs)
    //R -> CIC Filter Ratio
    //N -> CIC Filter Cascaded Sections
    //D -> Decimation Ratio
    w = w / D;
    h = abs(1/R*(sin(2*%pi*w*R/2) / sin(2*%pi*w/2)))^N;         //CIC Response
endfunction

function h = InverseCIC(w, R, N, D)
    //w -> Normalized frequncy point (1/Fs)
    //R -> CIC Filter Ratio
    //N -> CIC Filter Cascaded Sections
    //D -> Decimation Ratio
    w = w / D;
    if(w ~= 0) then
        h = abs(R*(sin(2*%pi*w*1/2) / sin(2*%pi*w*R/2)))^N;      //Inverse CIC Response
        //h = abs(M*(sin(2*%pi*1/(2*M)) / sin(2*%pi*w*1/2)));  //CIC Response
    else
        h = 1;
    end
endfunction

function h = LowPassWindow(w, f0, B)
    //w -> Normalized frequncy point (1/Fs)
    //f0 -> Cut-off frequency
    //B -> Transition bandwidth
    if(w <= f0) then
        h = 1;
    elseif(w <= (f0 + B)) then
        h = (1 - ((w - f0) / B));
    else
        h = 0;
    end
endfunction

function [c,w,hk] = FrequencySamplig(h, n, t)
    //h -> Desired frequency response function
    //n -> Number of coefficents to generate
    //t -> Sampling Type: either 1 or 2
    
    //Check if even or odd number of coeficients
    even = %t;
    if(modulo(n, 2) == 0) then
        even = %t;
    else
        even = %f;
    end
    
    //Get start and stop points for sampling symetry
    if(t == 1 && even == %t) then
        st = 1;
        sp = (n / 2) - 1;
    elseif(t == 1 && even == %f) then
        st = 1;
        sp = ((n - 1) / 2);
    elseif(t == 2 && even == %t) then
        st = 0;
        sp = (n / 2) - 1;
    elseif(t == 2 && even == %f) then
        st = 0;
        sp = ((n-1) / 2) - 1;
    end
    
    //Sample symmetric points
    for k = st:sp
        if(t == 1) then
            //Type 1 sampling
            wk = (2*%pi)/n * k;
        else
            //Type 2 sampling
            wk = %pi/n * (2*k + 1);
        end
        w(k+1) = exp(%i*wk);
        w((n-k)+st) = exp(-%i*wk);
        
        //Using complex sampling (aka with phase information)
        //ph = - modulo((2*%pi)/n * k, %pi);
        //mg = h(wk/(2*%pi));
        //hk(k+1) = mg*cos(ph) + %i*mg*sin(ph);
        //hk((n-k)+st) = mg*cos((-%pi - ph)) + %i*mg*sin((-%pi - ph));
        
        hk(k+1) = h(wk/(2*%pi));
        hk((n-k)+st) = h(wk/(2*%pi));
    end
    
    //Get missing points, the not symetric ones
    if(t == 1 && even == %t) then
        //H(0)
        k = 0;
        wk = (2*%pi)/n * k;
        w(k+1) = exp(%i*wk);
        
        //Using complex sampling (aka with phase information)
        //ph = - modulo((2*%pi)/n * k, %pi);
        //mg = h(wk/(2*%pi));
        //hk(k+1) = mg*cos(ph) + %i*mg*sin(ph);
        
        hk(k+1) = h(wk/(2*%pi));
        
        //H(N/2)
        k = n / 2;
        wk = (2*%pi)/n * k;
        w(k+1) = exp(%i*wk);
        
        //Using complex sampling (aka with phase information)
        //ph = - modulo((2*%pi)/n * k, %pi);
        //mg = h(wk/(2*%pi));
        //hk(k+1) = mg*cos(ph) + %i*mg*sin(ph);
        
        hk(k+1) = h(wk/(2*%pi));
    elseif(t == 1 && even == %f) then
        //H(0)
        k = 0;
        wk = (2*%pi)/n * k;
        w(k+1) = exp(%i*wk);
        
        //Using complex sampling (aka with phase information)
        //ph = - modulo((2*%pi)/n * k, %pi);
        //mg = h(wk/(2*%pi));
        //hk(k+1) = mg*cos(ph) + %i*mg*sin(ph);
        
        hk(k+1) = h(wk/(2*%pi));
    elseif(t == 2 && even == %t) then
        //None
    elseif(t == 2 && even == %f) then
        //H(N/2)
        k = (n - 1) / 2;
        wk = %pi/n * (2*k + 1);
        w(k+1) = exp(%i*wk);
        
        //Using complex sampling (aka with phase information)
        //ph = - modulo((2*%pi)/n * k, %pi);
        //mg = h(wk/(2*%pi));
        //hk(k+1) = mg*cos(ph) + %i*mg*sin(ph);
        
        hk(k+1) = h(wk/(2*%pi));
    end
    
    //Calculate coeficients
    for i = 0:(n-1)
        sumHk = 0;
        for k = 0:(n-1)
            wk = (2*%pi)/n * k;
            //sumHk = sumHk + (hk(k+1) * exp(%i*wk*i));
            sumHk = sumHk + real(hk(k+1) * exp(%i*wk*i));
        end
        c(i+1) = 1/n * sumHk;
    end

    if(even == %t) then
        c = circshift(c, floor(n/2) - 1);
        //c = [flipdim(c(1:ceil(n/2)), 1); c(1:ceil(n/2))];
    elseif(even == %f) then
        c = circshift(c, floor(n/2));
        //c = [flipdim(c(1:ceil(n/2)), 1); c(2:ceil(n/2))];
    end
endfunction

function c = CICCompFilter(m, d, fc, bp, n)
    //m -> CIC Filter Ratio
    //d -> Decimation Ratio
    //fc -> Low Pass cut-off frequency
    //n -> Number of coeficients/taps
    
    function h = CombinedFilter(w)
        h = InverseCIC(w, m, 1, d) * LowPassWindow(w, fc, bp);
    endfunction
    
    fr = linspace(0, 0.5, 512)';
    for i = 1:length(fr)
        xmCIC(i) = CIC(fr(i), m, 1, d);
        xmInvCIC(i) = InverseCIC(fr(i), m, 1, d);
        xmLP(i) = LowPassWindow(fr(i), fc, bp);
        xmComb(i) = CombinedFilter(fr(i));
    end
    
    //Mirror to get otherside
    xmCIC = [xmCIC; flipdim(xmCIC,1)];
    xmInvCIC = [xmInvCIC; flipdim(xmInvCIC,1)];
    xmLP = [xmLP; flipdim(xmLP,1)];
    xmComb = [xmComb; flipdim(xmComb,1)];
    
    //Actual frequncy sampling to design the compensation FIR filter
    [c, smp, hks] = FrequencySamplig(CombinedFilter, n, 1);
    
    //Plot graphs in Window
    scf(0);
    clf(0);
    
    //Plot CIC filter response
    fr = linspace(0, 1, 1024)';
    p = plot('nn', fr, xmCIC);
    p.foreground = color(128,0,255);
    //Plot inverse CIC filter response
    p = plot('nn', fr, xmInvCIC, ":");
    p.foreground = color(128,0,255);
    //Plot LPF window
    p = plot('nn', fr, xmLP, ":");
    p.foreground = color(255,255,255);
    //Plot ideal combined frequency response
    p = plot('nn', fr, xmComb, "g");
    p.foreground = color(0,255,0);
    
    //Plot actual designed filter response
    [xm, fr2] = frmag(c, 512);
    xm = [xm flipdim(xm,2)];
    p = plot('nn', fr, xm');
    p.foreground = color(255,0,0);
    
    //Draw sample points on top
    x = atan(imag(smp), real(smp)) / (2*%pi);
    for i = 1:length(x)
        if(x(i) < 0) then
            x(i) = 1 + x(i);
        end
    end
    p = plot('nn', x, hks, 'x');
    p.mark_foreground = color(255,0,0);
    
    //Draw Global Response, CIC with compensation
    xmNew = xm .* xmCIC';
    p = plot('nn', fr, xmNew');
    p.foreground = color(0,0,255);
    
    //Format plot3
    title("CIC Compensation Filter", 'fontname', "helvetica bold", 'fontsize', 6);
    a = get("current_axes");
    a.data_bounds = [0;0.5;0;1.2];
    a.x_ticks = tlist(["ticks" "locations", "labels"], [0.00; 0.05; 0.10; 0.15; 0.20; 0.25; 0.30; 0.35; 0.40; 0.45; 0.50], ['0.00'; '0.05'; '0.10'; '0.15'; '0.20'; '0.25'; '0.30'; '0.35'; '0.40'; '0.45'; '0.50']);
    a.y_ticks = tlist(["ticks" "locations", "labels"], [1.20; 1.00; 0.8; 0.60; 0.40; 0.20; 0.00], ['1.20'; '1.00'; '0.80'; '0.60'; '0.40'; '0.20'; '0.00']);
    a.sub_ticks = [0, 0];
    a.font_size = 4;
    a.font_style = 8;
    a.tight_limits = "on";
    xlabel("Frequency (1/Fs)", 'fontname', "helvetica bold", 'fontsize', 4);
    ylabel("Gain (x1)", 'fontname', "helvetica bold", 'fontsize', 4);
    xgrid();
    
    l = legend(['CIC';'Inverse CIC';'LPF Window';'Desired Response';'Compensation Filter';'FIR Samples';'Global Response'], 'in_upper_right');
    l.line_width = 0.05;
    l.background = color(0,0,0);
    l.foreground = color(255,255,255);
    l.font_color = color(255,255,255);
    l.font_size = 3;
    l.font_style = 8;
    
    //Syling
    hf = gcf();
    ha = gca();
    hf.background = color(0,0,0);           //Outer Area
    ha.background = color(0,0,0);           //Inner Area
    ha.foreground = color(255,255,255);     //Axes Rectangle
    ha.thickness = 4;                       //Line thickness
    ha.font_color = color(255,255,255);     //Ticks
    ha.grid = [color(255,255,255), color(255,255,255)];    //Grid
    ha.title.font_foreground = color(255,255,255);         //Title
    ha.x_label.font_foreground = color(255,255,255);       //x-Label
    ha.y_label.font_foreground = color(255,255,255);       //y-Label
    ha.children.children.thickness = 4;     //Line thickness
    
    xgrid(color(255,255,255), 1, 8);
endfunction
