(13 votes, average: 4.08 out of 5)

# Example of Cascaded Integrator Comb filter in Matlab

by on July 1, 2007

## 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 $N$ samples of an input sequence $x(n)$.

Figure: Moving average filter

The frequency response of the moving average filter is:

$H_f = \left|\frac{\sin(\pi fN)}{\sin(\pi f)}\right|$.

% 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 came$N$prior. 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:

$y(n) = \Sigma_{k=0}^{N-1}x(n-k) = x(n) - x(n-N) + y(n-1)$

#### 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 $D$ 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 $N$ which is achieved by repeating the sample $N$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

D id you like this article? Make sure that you do not miss a new article by subscribing to RSS feed OR subscribing to e-mail newsletter. Note: Subscribing via e-mail entitles you to download the free e-Book on BER of BPSK/QPSK/16QAM/16PSK in AWGN.

venkata sravani December 17, 2012 at 12:47 am

thank you very much for your explanation.

Hyeukjin,Lim July 29, 2012 at 11:18 am

delay buffer depth??
How do you decide ‘delay buffer depth’??

Krishna Sankar August 29, 2012 at 5:55 am

@Hyeukjin,Lim: Based on the filter cut off bandwidth required. Larger delay buffer depth translates to a lower cut off …

Prashant. S July 17, 2012 at 2:00 pm

Dear All,

I would like to know if there is some one who can show me an Matlab example which performs the Comb Filter Logic which is explained here in the below links. Any one logic will be enough to start with. This will help me to go forward with the filter implementation. It is very urgent for me as I have to meet my deadlines.

https://www.dropbox.com/s/5nhmxqgveytl0ba/Comb%20FIR%20Filter%20Explanation.jpg

(or)

https://www.dropbox.com/s/5nhmxqgveytl0ba/Comb%20FIR%20Filter%20Explanation.jpg

thanking you in advance for understanding!

Prashant

Prashant. S July 13, 2012 at 12:13 pm

Dear Krishna,

I want to know if I want to share files in this site then how can I do that? For example when I have a problem in MATLAB coding and I want to share file to show the people how I am coding it, as some times the codes are too long to paste here then it makes simpler to share the files instead.

Greetings

Krishna Sankar July 14, 2012 at 6:23 am

@Prashant: Well, you can upload to a public place like dropbox (or something equivalent) and share the link

Prashant.S July 15, 2012 at 12:59 am

Dear Krishna,

Here I have a FIR Comb Filter *.m file but I don’t know how I can modify it to meet my requirements as explained in my previous posts.

You need to just run TS.m file and you will see output of a Comb filter. I don’t if that will be helpful to both of my questions.

Greetings
Prashant

Krishna Sankar July 18, 2012 at 6:36 am

@Prashant: Knowing the sampling frequency of your systems, one can find where the nulls of the comb are present.

Prashant. S July 16, 2012 at 7:37 pm

Hallo Krishna,

I would like to know if there is a possibility to speak to you about my issue? If that is possible then it will be very helpful as I am getting close to my deadlines. Kindly let me know by what means will it be possible?

Sincerely
Prashant

Krishna Sankar July 18, 2012 at 6:37 am

@Prashant: Please look at my previous replies. Hope that helps

Prashant. S July 12, 2012 at 5:22 pm

Dear Krishna,

Apart from my previous post I have also the same Issue with a different requirement. Below are the link which shows what I am looking for. But the problem is how do I implement that? Here in this case I want to see only odd harmonics only in the output transfer function. I have never worked with ‘Comb Type FIR Filters’ so I am lost here.

I will be very grateful to you if you have a solution for such problem.
http://img821.imageshack.us/img821/7016/combfifofirfilter.jpg

The Below link shows how the buffer length ‘K’, Delay element Z^-1 & transmitted frequencies (Transfer Function) vary.

http://img826.imageshack.us/img826/4497/combfifofirfiltertabel.jpg

It will be very helpful in my project if you have any solution for such problem.

Thanking you!

Sincerely
Prashant.S

Krishna Sankar July 18, 2012 at 6:33 am

@Prashant: To get a comb response, one can use a circuit having a response as
y[n] = x[n] – x[n-D], where D is the delay

Putting D = 16, we can get a frequency response as follows:
B = [1 zeros(1,15) -1];
A = 1;
H W] = freqz([1 zeros(1,15) -1], [1],1024,’twosided’);
plot([-512:511]/1024,10*log10(abs(fftshift(H))));

Hope this helps.

Prashant. S July 18, 2012 at 12:35 pm

Hi krishna,

So I think ‘D’ will be order of the filter, right?
How can I mention the sampling freqency ‘Fs’ of this filter?

If your code is for this: https://www.dropbox.com/s/5nhmxqgveytl0ba/Comb%20FIR%20Filter%20Explanation.jpg

What part of this code tells me the following:
If Fs=1KHz, Input signal Period=1 then,
Buffer Length= (Input signal period/2)* Fs = (1/2)*Fs=500
There for:
First sample coming out of the Buffer=Buffer Length/Fs=500/1KHz=0.5 sec

Kindly guide me How I can determine these from the filter code.

Sincerely
Prashant

Krishna Sankar July 20, 2012 at 6:45 am

@Prashant: The concept of sampling frequency is notional in matlab. If we assume that sampling frequency is fs=1kHz, for the example with D=16, the nulls of the comb filter lies at +/-125Hz, +/-250Hz, +/-375Hz
B = [1 zeros(1,15) -1];A=1;
[H W] = freqz(B,A,1024,’twosided’);
plot([-512:511]/1024*1000,10*log10(abs(fftshift(H))));
xlabel(‘freq, Hz’);ylabel(‘amplitude, dB’);grid on
octave:22> axis([-500 500 -15 5])

liguo May 31, 2012 at 11:01 am

i want to plot the magnitude response of CIC filter.Thanks in advance.That would help me a lot.

Krishna Sankar July 2, 2012 at 5:18 am

@liguo: You can plot the frequency response of the moving average filter as
% 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’)

arunraj April 26, 2012 at 7:33 pm

sir i need exact defination of cic filter……need of cic flter ,working ,advantages and disadvantages of cic filter……………..please send it as soon as possible

Krishna Sankar April 27, 2012 at 5:46 am

@arunraj: does this post help?

Bhavani December 22, 2011 at 12:01 pm

i want to plot the magnitude response of comb filter…..anyone help me with Matlab code…Thanks in advance………cic fiter

Krishna Sankar January 4, 2012 at 5:37 am

@Bhavani: Hope this matlab code snippet helps.
% freq response
h_comb = [1 zeros(1,7) -1];
hF = fft(h_comb,1024);
plot([-512:511]/1024, 10*log10(abs(fftshift(hF))))
xlabel(‘normalized frequency’); ylabel(‘amplitude, dB’); title(‘Amplitude response of comb filter’);
grid on;

Alexander April 27, 2011 at 6:53 pm

I’m very interested about the next:
You said:
“Considering so, do we need to have the filtering hardware when doing interpolation? I do not think so. ”

Is there a way to replace the block “arrow up, I”,with zero-hold block(which passes the value of the input signal to the output instead of inserting a zero between samples).

How i can realize that in matlab?

Krishna Sankar May 26, 2011 at 5:43 am

@Alexander: You can insert zeros and then convolve with a rectangular window
For eg,
xt = [1 zeros(1,2) 2 zeros(1,2) 3 zeros(1,2)];
ht = ones(1,3);
yt = conv(xt,ht);

richa November 9, 2010 at 6:00 pm

hello krishna
nice to see your article on cic decimation filter
i am having problem with frequency response of cic specifically i want to improve it
can you help me i will be highly obiliged

Krishna Sankar November 14, 2010 at 10:34 am

@richa: Please provide more details of the problem which you are facing.

jennifer March 16, 2010 at 2:33 pm

hello..
We are doing our project on uwb signal in nano seconds ..we need code for down sampling that received wideband signal after bandpass filtering.plz help us sir

Krishna Sankar March 28, 2010 at 2:32 pm

@jennifer: Sorry, I have not tried modeling UWB signals. However,
for downsampling by factor of 2, take one sample every two samples;
for downsampling by factor of 3, take one sample every three samples; and so on

mike_cv February 7, 2010 at 10:19 am

hi. i had a homework to do but i can’t do it. so please i need your help to solve it. it’s about matlab and digital signal processing.

Krishna Sankar April 4, 2010 at 4:16 am

@mike_cv: What is the help which you are looking for

thomas August 25, 2009 at 10:56 am

Hi Krishna,

I still can not understand whether we must use digital low pass filter after interpolation. In fact , for upsampling, there is no aliasing issue at all and also the spectrum of the interpolation sequence is still periodic through low pass filter. To my understanding, after interpolation without low pass filter we can deliver the signal to DAC module directly to get the analog signal in time domain.Right ?

In addition, is there any other way to understand whether we can move the Decimation from the end of the cic to the mid of the cic without the help of matlab? That is to say, could we understand this by formula?

Krishna Sankar September 7, 2009 at 4:52 am

@thomas: My replies:
1/ Any digital signal will have spectrum from -infinity to +infinity. So, if we do not use the low pass filter, we will be transmitting the original spectrum plus all the replicas. And this is not desirable.
Yes, one can drive the output of the interpolation (inserting of zeros between samples) and pass to the DAC. However, note that most DAC’s have a zero-order hold transfer function i.e. hold on to the previous sample till the next sample has arrived. This zero-order hold is simple rectangular filter with sinc shaped frequency response.

2/ I have not tried playing with the math. However, it seems that once we write the transfer function of each stage, we should be able to get the equivalence.

Julián July 5, 2009 at 7:17 pm

Hi!

One question. If I have to create harmonics of a audio signal… I want make a basic integrator, that it have harmonic distorsion. Can you help?
I want a filter coefficients (A and B)for a Matlab function filter(B,A,x)!!!
Thanks!
PS: Sorry for my English.

Krishna Sankar July 6, 2009 at 7:21 pm

@Julián: A simple integrator can be a filter with transfer function H(z) = kz^-1/(1 – (1-k)z^-1)
% example matlab code
ip = ones(1,100);
k = 0.2;
op2 = filter([0 k],[1 -(1-k)],ip);

samanvita March 25, 2009 at 1:27 pm

Hi,
could you please explain the frequency response code snippet given above in one of the posts. How does it differ for a decimation filter?

regards,
Samanvita

Krishna Sankar April 4, 2009 at 7:31 am

@samanvita: Whether its interpolation or decimation filter, both are rectangular filters in time domain, which has a sinc() shaped frequency response. The rectangular filters can be equivalently implemented using an integrator arm and a comb arm . In the case of differentiation filter, we have the comb arm after the downsampling (to save the delay line depth).

Hope this helps.

Rohith November 18, 2008 at 12:10 am

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

VA April 12, 2010 at 9:40 am

hello,
In this post krishna has worked on single stage Integrator and comb.
and its true that in Hogenauer’s paper – “An economical class of digital filters for decimation and interpolation” he said about multistage implementation.
As we increase the number of stages the image/alias rejection will improve.
Here in krishna blog N=10 only refers to the number of samples in T duration.
I hope this will help.

Krishna Sankar November 15, 2008 at 8:12 am

@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.

Rohith November 12, 2008 at 6:18 am

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.”.

Krishna Sankar November 9, 2008 at 8:01 pm

@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.

Rohith November 8, 2008 at 1:50 am

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

Gary October 30, 2008 at 12:11 am

This is great! Thanks!

Krishna Sankar April 30, 2008 at 7:37 am

@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.

Rudheesh Raghav April 30, 2008 at 6:57 am

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??

Krishna Sankar March 4, 2008 at 5:51 am

@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.

venkata March 4, 2008 at 5:21 am

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

Prashant. S July 12, 2012 at 4:49 pm

Dear Venkat,

Though it is long time ago you have posted. I just want to give you the code for fixed point FIR Lowpass filter in Matlab:

Fs = 1400; % Sampling Frequency
Fc =400; % Cutoff Frequency
N = 20; % Order of the Filter
n = 8; % Input word length i.e if ADC is 8 bit then use 8, if 12 bit then use %12,…and so on
CWL = 8; % Coeff word length must be >= to ‘n’ so i.e if n is 8 bit then use 8 or %>8, if 12 bit then use 12 or >12,…and so on

flag = ‘scale’; % Sampling Flag
% Create the window vector for the design algorithm.
win = rectwin(N+1);
% Calculate the coefficients using the FIR1 function.
b = fir1(N, Fc/(Fs/2), ‘low’, win, flag);
Hd = dfilt.dffir(b);
% Set the arithmetic property.
set(Hd, ‘Arithmetic’, ‘fixed’, …
‘CoeffWordLength’, CWL, …
‘CoeffAutoScale’, true, …
‘Signed’, true, …
‘InputWordLength’, n, …
‘inputFracLength’, 0, …
‘FilterInternals’, ‘FullPrecision’);
denormalize(Hd);

% plots the magnitude response of the filter by default. If you have %a signal %processing tool box then you can explore more on this plot

fvtool(Hd)

%**********************************************************************%
%***Below code is to convert double type coefficients to integer****%
%******Write filter coefficients to GUI Table/excel sheet**************%
%*********************************************************************%
conv_coefficients_MAT=cell2mat(conv_coefficients_CELL);
Coeff_MAX_Val=max(abs(conv_coefficients_MAT));
Integer_Coeff_Val=round(((2^(CWL-1)-1)*conv_coefficients_MAT/Coeff_MAX_Val));
xlswrite(‘Lowpass_Filter_Coefficients.xls’, (Integer_Coeff_Val)’);

Hope that will help.

Prashant.S

Krishna Sankar July 13, 2012 at 5:23 am

@Prashant: Thanks for the help.

Krishna Sankar February 21, 2008 at 7:54 am

@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’)

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

TULASI BABU D February 20, 2008 at 3:51 pm

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

Bulletin News December 19, 2007 at 1:35 am

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

Krishna Sankar October 31, 2007 at 8:06 am

@ 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 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

priya October 30, 2007 at 4:57 pm

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

Krishna August 8, 2007 at 10:12 am

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

Anonymous August 8, 2007 at 1:18 am

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.