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
.
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’)
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:
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.
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).
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.
If you liked this post, you may leave a comment below, or subscribe to the RSS feed.
You may also find these posts relevant...
Comments
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?
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
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??
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
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.”.
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

































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.