***17

Example of Cascaded Integrator Comb filter in Matlab

Equivalence of Moving Average and CIC filter

Let me briefly share my understanding on the cascaded integrator comb (CIC) filter, thanks to the nice article. For understanding the cascaded integrator comb (CIC) filter, firstly let us understand the moving average filter, which is accumulation latest samples of an input sequence .


Moving Average filter

Figure: Moving average filter

The frequency response of the moving average filter is:

.

% Moving Average filter
N = 10;
xn = sin(2*pi*[0:.1:10]);
hn = ones(1,N);
y1n = conv(xn,hn);

% transfer function of Moving Average filter
hF = fft(hn,1024);
plot([-512:511]/1024, abs(fftshift(hF)));
xlabel(’Normalized frequency’)
ylabel(’Amplitude’)
title(’frequency response of Moving average filter’)

cic filter frequency response

Figure: Frequency response of moving average filter

The moving average filter which is implemented as a direct form FIR type as shown above can also be implemented in a recursive form. It consists of a comb stage whose output is difference of the current sample and the sample which cameprior. The difference is successively accumulated by an integrator stage. Together the circuits behave identically as the moving average filter described prior.

It can be proved using a simple proof that:

cascaded integrator comb filter

Figure: Cascaded Integrator Comb Filter

Bit more details can be found in Sec2.5.2 Recursive and Non Recursive Realization of FIR systems in [DSP-PROAKIS].

As the system is linear, the position of the integrator and comb stage can be swapped. Now, using a small Matlab code snippet let us verify that the output from CIC realization is indeed the same as obtained from moving average filter.

% Implementing Cascaded Integrator Comb filter with the
% comb section following the integrator stage
N = 10;
delayBuffer = zeros(1,N);
intOut = 0;
xn = sin(2*pi*[0:.1:10]);
for ii = 1:length(xn)
% comb section
combOut = xn(ii) - delayBuffer(end);
delayBuffer(2:end) = delayBuffer(1:end-1);
delayBuffer(1) = xn(ii);

% integrator
intOut = intOut + combOut;
y2n(ii) = intOut;
end

err12 = y1n(1:length(xn)) - y2n;
err12dB = 10*log10(err12*err12′/length(err12)) % identical outputs
close all

% Implementing Cascaded Integrator Comb filter with the
% integrator section following the comb stage

N = 10;
delayBuffer = zeros(1,N);
intOut = 0;
xn = sin(2*pi*[0:.1:10]);
for ii = 1:length(xn)
% integrator
intOut = intOut + xn(ii);

% comb section
combOut = intOut - delayBuffer(end);
delayBuffer(2:end) = delayBuffer(1:end-1);
delayBuffer(1) = intOut;
y3n(ii) = combOut;

end
err13 = y1n(1:length(xn)) - y3n;
err13dB = 10*log10(err13*err13′/length(err13)) % identical outputs

The outputs are matching.

The recursive realization of the FIR filter as described above helps to achieve the same result with less hardware.

Using CIC filter for decimation

Typically, decimation to a lower sampling rate is achieved by taking one sample out of every samples. There exists an anti-aliasing filter to remove the un-desired spectrum before decimation.

cic filer for decimation

Figure: CIC filters for decimation

As shown above, the same output can be achieved by having the decimation stage between integrator stage and comb stage. This helps in reducing the delay buffer depth requirement of the comb section. Using a small Matlab code snippet, let us check whether both the implementations behave identically.

% For decimation, having the CIC filtering before taking every other sample
D = 2; % decimation factor
N = 10; % delay buffer depth
delayBuffer = zeros(1,N); % init
intOut = 0;
xn = sin(2*pi*[0:.1:10]);
y6n = [];
for ii = 1:length(xn)
% comb section
combOut = xn(ii) - delayBuffer(end);
delayBuffer(2:end) = delayBuffer(1:end-1);
delayBuffer(1) = xn(ii);

% integrator
intOut = intOut + combOut;
y6n = [y6n intOut];

end
y6n = y6n(1:D:end); % taking every other sample - decimation

% For efficient hardware implementation of the CIC filter, having the
% integrator section first, decimate, then the comb stage
% Gain : Reduced the delay buffer depth of comb section from N to N/D
D = 2; % decimation factor
N = 10; % delay buffer depth
delayBuffer = zeros(1,N/D);
intOut = 0;
xn = sin(2*pi*[0:.1:10]); % input
y7n = []; % output
for ii = 1:length(xn)
% integrator
intOut = intOut + xn(ii);

if mod(ii,2)==1
% comb section
combOut = intOut - delayBuffer(end);
delayBuffer(2:end) = delayBuffer(1:end-1);
delayBuffer(1) = intOut;
y7n = [ y7n combOut];
end

end
err67 = y6n - y7n;
err67dB = 10*log10(err67*err67′/length(err67))

The outputs are matching.

Using CIC filters for interpolation

Typically, interpolation to a higher sampling rate achieved by inserting zeros between consecutive samples followed by filtering (for removing the images).

cic filter for interpolation

Figure: Using CIC filters for interpolation

As shown above, the same result can be achieved by having the upsampling stage between comb and integrator stage. This helps in reducing the delay buffer depth requirement of the comb section. Using a small Matlab code snippet, let us check whether both the implementations behave identically.

% For interpolation, insert the zeros followed by CIC filtering
xn = sin(2*pi*[0:.1:10]);
I = 2; % interpolation factor
N = 10; % filter buffer depth
xUn = [xn; zeros(1,length(xn))];
xUn = xUn(:).’; % zeros inserted
delayBuffer = zeros(1,N);
intOut = 0;
for ii = 1:length(xUn)
% comb section
combOut = xUn(ii) - delayBuffer(end);
delayBuffer(2:end) = delayBuffer(1:end-1);
delayBuffer(1) = xUn(ii);

% integrator
intOut = intOut + combOut;
y4n(ii) = intOut;

end

% For efficient hardware implementation of CIC filter for interpolation, having
% the comb section, then zeros insertion, followed by integrator section
% Gain : Reduced the delay buffer depth of comb section from N to N/I

I = 2; % interpolation factor
N = 10; % original delay buffer depth
delayBuffer = zeros(1,N/I); % new delay buffer of N/I
intOut = 0;
xn = sin(2*pi*[0:.1:10]);
y5n = [];
for ii = 1:length(xn)
% comb section
combOut = xn(ii) - delayBuffer(end);
delayBuffer(2:end) = delayBuffer(1:end-1);
delayBuffer(1) = xn(ii);

% upsampling
combOutU = [ combOut zeros(1,1)];

for jj =0:I-1
% integrator
intOut = intOut + combOutU(jj+1);
y5n = [y5n intOut];
end

end

err45 = y4n - y5n;
err45dB = 10*log10(err45*err45′/length(err45)) % outputs matching

The outputs are matching.

One one question: Upsampling by a factor which is achieved by repeating the sample times result in the same output as obtained by the CIC filter implementation described above (see previous post trying to describe that zero-order hold for interpolation). Considering so, do we need to have the filtering hardware when doing interpolation? I do not think so. If I have some additional thoughts, I will update.

References:

[1] Understanding cascaded integrator-comb filters - By Richard Lyons, Courtesy of Embedded Systems Programming URL: http://www.us.design-reuse.com/articles/article10028.html

[2] Digital Signal Processing - Principles, Algorithms and Applications, John G. Proakis, Dimitris G. Manolakis

Please click here to SUBSCRIBE to newsletter and download the FREE e-Book on probability of error in AWGN. Thanks for visiting! Happy learning.

1 Star2 Stars3 Stars4 Stars5 Stars (No Ratings Yet)
Loading ... Loading ...

If you liked this post, you may leave a comment below, or subscribe to the RSS feed.

You may also find these posts relevant...
  • Articles
  • Sigma delta modulation
  • Update: Correction in Matlab code for raised cosine filter
  • Interpreting the output of fft() operation in Matlab
  • Raised cosine filter for transmit pulse shaping
  • Comments

    For moving average,
    |Hf|
    ={1/N*[sin(pi*f*N)]/[sin(pi*f)]}

    but for CIC
    |Hf|
    ={[sin(pi*f*N)]/[sin(pi*f)]}

    Why they are the same, as a scaling factor is there, 1/N ? thanks.

    Yes, for notational simplicity I ignored the scaling factor 1/N. However, I do believe that you agree with the recursive implementation.
    With the scaling factor,
    y[n] = y[n-1] + (1/N)*(x[n] - x[n-N])

    % simple matlab example:
    ip = randn(1,100);
    op1 = conv(ip,(1/16)*ones(1,16));
    op2 = filter((1/16)*[1 zeros(1,15) -1],[1 -1],ip);
    err = op1(1:100) - op2;
    errdB = 10*log10(err*err’/length(err))

    HTH,
    Krishna

    how to plot cic filter frequency response in matlab with the code mentioned above for interpolation?

    @ priya

    Background:
    As you may have observed, the CIC filter is a recursive implementation of a moving average filter (let us say of length N). Typically, when interpolation by a factor I , we insert I-1 zeros inbetween the samples and then do a filtering to remove the aliases.

    When doing interpolation, for efficient hardware implementation (i.e. to reduce the number of elements in the comb buffer), we can first do the comb (with N/I delay elements), insert I-1 zeros, then integrate. However, the effect is same as passing the I-1 zero inserted samples through a moving average filter of length N.

    To answer your query:
    To find the frequency response, you can either use fft() of the time domain impulse response, or use the freqz() function in Matlab. Both gives identical frequency response.

    % simple example - matlab/octave
    clear all
    close all
    N = 16;
    fs = 1;
    ht = ones(1,N);
    h1F = fft(ht,1024);
    [h2F f2] = freqz([1 zeros(1,N-1) -1],[1 -1],1024,’whole’);
    plot([-512:511]*fs/1024,abs(fftshift(h1F)),’b.-’)
    hold on
    plot([-512:511]*fs/1024,abs(fftshift(h2F)),’rx-’)
    xlabel(’frequency’)
    ylabel(’magnitude’)

    Did that help?

    Krishna

    Amazing blog post about le of Cascaded Integrator Comb filter in Matlab - DSP log. I enjoy this view!

    dear sir,
    this is tulasi
    am working on “comb generator” i know it basically under microwave product.

    my basic doubt is 1) what is comb generator?

    2) what is the difference between comb generator and comb filter?
    3) i need matlab implementation of comb generator software code

    thanking u
    tulasi babu

    hyderabad

    @tulasi:
    the ‘comb’ aspect in CIC filter corresponds to the frequency response of the ‘comb’ section. If you use matlab/octave and run the attached code snippet, you can see that frequency response of the section
    y[n] = x[n] - x[n-16] looks like start-quote “The periodic zeros remind us of the teeth of a comb, hence the name, comb filter.” end-quote

    clear all
    close all
    N = 16;
    fs = 1;
    [h1F f1] = freqz([1 zeros(1,N-1) -1],[1],1024,’whole’);
    plot([-512:511]*fs/1024,abs(fftshift(h1F)),’b.-’)
    hold on
    xlabel(’frequency’)
    ylabel(’magnitude’); title(’frequency response of comb section’)

    Hope this answers your query#2.

    For query#1, from a quick wiking, I found that comb generator is a source which can generate multiple harmonics of the same signal.

    and query#3: generating waveform which is the sum of multiple harmonics of the same sinusoidal should be reasonably trivial, no? You can look at the OFDM posts.

    HTH,
    Krishna
    ps. it is not a good idea to provide mobile phone numbers on public domain. hence i edited that out from the post

    can anyone got the code for designing theFIR lowpass filter in MATLAB

    @venkata : you can use the fdatool for designing both fir/iir low/high/band pass/stop filters in Matlab.
    type fdatool on the matlab command line to invoke the gui.

    Hi Krishna,

    Thank you for the excellent article. You have simplified the concept in a way so that evryone can understand.

    Keep up the good work.

    By the way, are you from kerala??

    @Rudheesh: Thanks. You can also also look at the other posts and give your valuable feedback.

    Yeah, am from Kerala…. though now its Bangalore which gives me the rice and sambar. ;)

    This is great! Thanks!

    Hello Krishna,

    I’m not sure if I completely follow the code for the CIC filter.

    Does “N” refer to the number of stages? If so, then why is it that the integrator only has one stage?

    Here is the code from the above post for CIC filter where the comb section follows the integrator section:

    Thanks,
    Rohith
    % Implementing Cascaded Integrator Comb filter with the
    % integrator section following the comb stage

    N = 10;
    delayBuffer = zeros(1,N);
    intOut = 0;
    xn = sin(2*pi*[0:.1:10]);
    for ii = 1:length(xn)
    % integrator
    intOut = intOut + xn(ii);

    % comb section
    combOut = intOut - delayBuffer(end);
    delayBuffer(2:end) = delayBuffer(1:end-1);
    delayBuffer(1) = intOut;
    y3n(ii) = combOut;

    end

    @Rohith: Its the comb stage which has the N delay element, integrator has only one delay. I am hoping that the figure above is self explanatory.

    Hello Krishna,

    Thanks for the reply. I noticed that in the diagram too, which is where I got confused. The code and the diagram go together, but I am wondering whether the concept is right.

    Here is a excerpt from Hogenauer’s original paper on CIC filters. Note that it mentions “N” stages for both integrator and comb sections.

    “The integrator section of CIC filters consists of N ideal digital integrator stages operating at the high sampling rate, fs,. Each stage is implemented as a one-pole filter with a unity feedback coefficient. The comb section operates at the low sampling rate fs/R, where R is the integer rate change factor. This section consists
    of N comb stages with a differential delay of M samples per
    stage.”.

    @Rohit: I downloaded the Hogenauer’s paper - “An economical class of digital filters for decimation and interpolation”

    As I understand, he is just stating there can be N integrator stages with each having unity delay. See Eq (1) and Fig (1,2).

    In our simple example, we are having only ONE integrator stage. Hope this explains.

    Hello Krishna,

    Thanks for looking into the paper. I just wanted to point out that “N” has a major impact on the filter. As you have noticed, Eq(1) refers to a single integrator stage and Eq(2) refers to a single comb stage. The number of stages in the integrator and the number of stages in the comb section are equal. The only difference is that the comb section runs at a lower sampling rate.

    Anyways, glad you were able to look into this paper. May be it’ll give you ideas for the future.

    Regards,
    Rohith

    Leave a comment

    (required)

    (required)