% CIC decimator example
% includes 
% - time domain implementation
% - frequency domain model via z-domain transfer function
% Parameters:
% R = rate change factor
% N = nStages
% M = differential delay
% reference: [1] http://www.altera.com/literature/an/an455.pdf 

function instantCIC()
    close all;
    testvec = zeros(1, 43);
    R = 4;
    testvec(1) = 1;
    N = 3;
    M = 2; 
    a = CICdec_timeDomainModel(testvec, N, M, R);
    [b, H] = CICdec_freqDomainModel(testvec, N, M, R, false);
    [c, H] = CICdec_freqDomainModel(testvec, N, M, R, true);

    figure(1); grid on; hold on;
    plot(1:R:R*numel(a), a, 'bx');
    plot(b, 'ro');
    plot(1:R:R*numel(c), c, 'k+');
    title('CIC decimator impulse response');
    legend('time domain', 'freq. domain', 'freq. domain decimated');

    % plot
    H = fft([ifft(H), zeros(1, 1000)]); % zero-padding
    figure(2); clf(); grid on; hold on; 
    plot(linspace(-0.5, 0.5, numel(H)), fftshift(20*log10(abs(H) + 1e-15)), 'b'); 
    xlim([-0.5, 0.5]);
    title('frequency response');
    xlabel('f/fs');
    ylabel('dB');
end

function [vec, H] = CICdec_freqDomainModel(vec, N, M, R, doDecim) 
    flag = isreal(vec); 
    
    % evaluate frequency response (z-domain transfer function)
    n = numel(vec);
    zInv = exp(-2i*pi*(0:(n-1))/n); 
    b = ones(1, R * M);
    H = polyval(b, zInv) .^ N;
    H = H / H(1);

    % apply frequency response
    vec = ifft(fft(vec) .* H);
    
    % decimate
    if doDecim
        vec = vec(1:R:end);
    end
    
    % don't let FFT roundoff error turn real signal into complex
    if flag
        vec = real(vec);
    end 
end

function vec = CICdec_timeDomainModel(vec, N, M, R)
    nLeadIn = M * N * R;
    nLeadOut = nLeadIn / R;
    ix = mod(-nLeadIn:-1, numel(vec)) + 1;
    
    % prepend end of cyclic signal
    vec = [vec(ix) vec];
    
    % integrator
    for ix = 1:N
        vec = cumsum(vec);
    end

    % decimator
    vec = vec(1:R:end);
    
    % differentiator
    for ix = 1:N
        vec = vec - circshift(vec, [0, M]); 
    end

    % remove the added length
    vec = vec(nLeadOut+1:end);

    % scale gain
    gain = (R * M) .^ N;
    vec = vec / gain;
end
