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

thank you very much for your explanation.

Now i understood completely about CIC filter and its advantage.

delay buffer depth??

How do you decide ‘delay buffer depth’??

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

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

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

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

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.

https://dl.dropbox.com/u/91710211/Comb%20Filter%20m%20Files/comb.m

https://dl.dropbox.com/u/91710211/Comb%20Filter%20m%20Files/TS.m

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

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

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

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

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

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

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

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

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

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

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

@arunraj: does this post help?

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

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

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?

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

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

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

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

@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

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.

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

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?

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

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.

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

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

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

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

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.

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

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

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

This is great! Thanks!

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

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

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

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

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**************%

%*********************************************************************%

read_coefficients = Hd.coeffs; %read_coefficients is structure type

conv_coefficients_CELL=struct2cell(read_coefficients);

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

@Prashant: Thanks for the help.

@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

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

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

@ 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

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

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

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.