Forums

Cascaded Integrator Comb and aliasing experts please help!

Started by David Joseph Bonnici June 5, 2004
Thanks for not losing your patience.  I am going to make a resume of
what my intentions are and list my two questions in order of priority.
 I receive your postings with a latency of 5 hours so if someone has
already made postings giving me solutions please discard these
message.

Question 1: I need to compute the aggregate frequency response of
three cascaded multirate filters.
I was trying the following but I was in mistake:
What I was trying to do was to get the impulse response of all the
multirate filters and the convolve them. (as suggested on the DDC
datasheet of Xilinx)  The problem was that I was not considering the
filters as multirate and so they were being referred to just one
sampling frequency (Fs) and this is not true since Fs changes.
After I have recognized that I had a problem I tried a different
solution:
I've inputted an impulse to just the CIC.  After getting the response
(which is also decimated) I've inputted the output waveform into the
PFIR and so on until I get the output from the CFIR.  The result I
have obtained does not match with published results of the GSM DDC
from Xilinx.  I have tested the all on Matlab R13 which does not have
support of the multirate functions.  I have used the following code.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
clc; %GSM Spectral Mask test so as to check the functionality of the
script with that published in Xilinx DDC
clear;
input = zeros(1000,1); input(1) = 1;   % An impulse

%%%%%%%%%%%%%%%%CIC%%%%%%%%%%%%%%%%%%%%%%%%%%%%
q = quantizer([8 0],'fixed'); % Defines the properties of input.
output_cic = cicdecimate(1, 4, 48, input, q);

figure;
plot(output_cic); title('Time plot for output due to impluse response
of the CIC');
freq_y = fft(output_cic,4096);
figure;
l=length(freq_y);
fr=0:1/(l-1)*100e6:100e6;
plot(fr,abs(freq_y)); xlabel('Frequency,Hz'); ylabel('Amplitude');
title('Frequency Response');
figure;
plot(fr,20*log10(abs(freq_y))); xlabel('Frequency,Hz');
ylabel('Amplitude in dB'); title('Frequency Response');
figure; plot(fr, unwrap(angle(freq_y))); xlabel('Frequency,Hz');
ylabel('Phase'); title('Phase for y');
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%PFIR%%%%%%%%%%%%%%%%%%%%%%%%
h = [-98,-679,-2016,-3234,-2537,850,6053,12060,18230,23239,25212,
23239,18230,12060,6053,850,-2537,-3234,-2016,-679,-98];

p =1;
q =2;
output_c = upfirdn(output_cic,h,p,q)
figure;
plot(output_c); title('C section');
freq_y = fft(output_c,4096);
figure;
l=length(freq_y);
fr=0:1/(l-1)*100e6:100e6;
plot(fr,abs(freq_y)); xlabel('Frequency,Hz'); ylabel('Amplitude');
title('Frequency Response');
figure;
plot(fr,20*log10(abs(freq_y))); xlabel('Frequency,Hz');
ylabel('Amplitude in dB'); title('Frequency Response');
figure; plot(fr, unwrap(angle(freq_y))); xlabel('Frequency,Hz');
ylabel('Phase'); title('Phase for y');
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


%%%%%%%%%%%%%CFIR%%%%%%%%%%%%%%%%%%%%%%
h = [1007,-1853,79,807,1633,423,265,175,-527,-1331,-1454,-1087,-721,-149,1008,1985,
2164,2005,1483,61,-1756,-3134,-3953,-4016,-2714,134,4003,8361,12934,17009,19565,20353,
19565,17009,12934,8361,4003,134,-2714,-4016,-3953,-3134,-1756,61,1483,2005,2164,1985,1008,-149,
-721,-1087,-1454,-1331,-527,175,265,423,1633,1807,79,-1853,1007];

p =1;
q =2;
output_p = upfirdn(output_c,h,p,q)
figure;
plot(output_p); title('P section');
freq_y = fft(output_p,4096);
figure;
l=length(freq_y);
fr=0:1/(l-1)*100e6:100e6;
plot(fr,abs(freq_y)); xlabel('Frequency,Hz'); ylabel('Amplitude');
title('Frequency Response');
figure;
plot(fr,20*log10(abs(freq_y))-293.1); xlabel('Frequency,Hz');
ylabel('Amplitude in dB'); title('Frequency Response');
figure; plot(fr, unwrap(angle(freq_y))); xlabel('Frequency,Hz');
ylabel('Phase'); title('Phase for y');
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

You might note that I have been using upfirdn so as to emulate the
polyphase filtering.  I do not have a better function so as to emulate
the polyphase filtering.  I do not have the latest version of Matlab
so I do not have access to better functions.  Still it is not supposed
to make difference since the upsampling ratio is 1 and I am supplying
the filter coefficients myself.

Now is the problem related with the use of upfirdn or is the
underlying taught completely wrong? What is the correct way of doing
it?

Question Number 2:
In order of importance Question 1 is more important than this one but
if you can highlight a possible solution I would appreciate it.
Andraka has already suggested a solution but I did not understand what
is he aiming to.
I have used the thing of inserting the zeros between the coefficients
for a different purpose.  I have used it so as to get the impulse
response of the CIC referenced to the high sampling frequency.
The code is the following
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clc; 
clear;
znum = zeros(193,1);
znum(1) = 1;
znum(49) = -4;
znum(97) = 6;
znum(145) = -4;
znum(193) = 1; 
zden = zeros(193,1);
zden(1) = 1;
zden(2) = -4;
zden(3) = 6;
zden(4) = -4;
zden(5) = 1;

input = zeros(200,1); input(1) = 1;   % An impulse
output_cic = filter(znum,zden,input); %the CIC is treated as an IIR

figure;
plot(output_cic); title('Time plot for output due to impluse response
of the CIC');
freq_y = fft(output_cic,4096);
figure;
l=length(freq_y);
fr=0:1/(l-1)*100e6:100e6;
plot(fr,abs(freq_y)); xlabel('Frequency,Hz'); ylabel('Amplitude');
title('Frequency Response');
figure;
plot(fr,20*log10(abs(freq_y))); xlabel('Frequency,Hz');
ylabel('Amplitude in dB'); title('Frequency Response');
figure; plot(fr, unwrap(angle(freq_y))); xlabel('Frequency,Hz');
ylabel('Phase'); title('Phase for y');
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
The code does its job correctly.

That's all.  Hope you can help me.
Many Thanks again
David Joseph Bonnici
Any ideas, because I am getting crazy about this.  I cannot believe
it, that there is no where cited how to get the composite frequency
response of cascaded mulitrate filters.  Is it,  because it is either
too simple or either too complicated?
It is practically impossible to design cascaded multirate filters
without seeing the aggregate frequency response.  I have observed a
too high degree risk of error when designing cascaded multirates.

David
David Joseph Bonnici wrote:

> Any ideas, because I am getting crazy about this. I cannot believe > it, that there is no where cited how to get the composite frequency > response of cascaded mulitrate filters. Is it, because it is either > too simple or either too complicated? > It is practically impossible to design cascaded multirate filters > without seeing the aggregate frequency response. I have observed a > too high degree risk of error when designing cascaded multirates. > > David
When all else fails, multiply the frequency responses (add the dB responses) frequency by frequency. It's relatively easy with analytic expressions, and only somewhat tedious if all you have is curves. Jerry -- Engineering is the art of making what you want from things you can get. �����������������������������������������������������������������������

David Joseph Bonnici wrote:
> > Any ideas, because I am getting crazy about this. I cannot believe > it, that there is no where cited how to get the composite frequency > response of cascaded mulitrate filters. Is it, because it is either > too simple or either too complicated?
Wasn't that what the ^N was all about in your original question? To get the frequency response of a series of cascaded filters you just multiply the frquency response of each. If they are all identical filters then its just a matter of raising the response to the power of N (N is the number of filters). -jim -----= Posted via Newsfeeds.Com, Uncensored Usenet News =----- http://www.newsfeeds.com - The #1 Newsgroup Service in the World! -----== Over 100,000 Newsgroups - 19 Different Servers! =-----
First, the CIC is actually an FIR filter expressed in a recursive form.  If the CIC has a
decimation ratio of R and a delay of N samples in the comb section, then the response is the same
as an N*R tap boxcar (all coefficients = 1) FIR filter referred to the higher sample rate.  In many
instances the CIC will be the first decimation stage followed by an FIR filter, or perhaps a series
of decimating FIR filters.  The CIC's response in the frequency domain is the fourier transform of
a step (because of the boxcar shape of the equivalent FIR filter), which is the sinc function
cic_resp=abs((sin(pi*n*r*f)./(n*r*sin(pi*f))).^m);  If you have a higher order CIC, then raise the
response to the power m.
The first FIR filter is running at a decimated rate.  If for example the CIC decimates by 8, then
the FIR sample rate is 1/8 of the CIC rate.  The response of an FIR filter in the frequency domain
is the fourier transform of the coefficients (the coefficients of an FIR are the impulse
response).  If you insert zero's between each FIR filter coefficient, then you get the original
frequency domain response but compressed in frequency and replicated (see what happens when you
insert 7 zeros between each coefficient and take the fourier transform of the new filter, compare
it to the original).  Note that this if you also expand the frequency axis markings by multiplying
by 8 you get frequency axis markings matching your original filter on the first copy of the
original response.  This expanded filter is therefore referred to the input sample rate at the CIC
input.  To get the composite response, use the log of the individual frequency responses and add
them together.

Here is a simple example:
%CIC response mapped to input side of CIC:
cic_resp=20*log10(abs((sin(pi*n*r*f)./(n*r*sin(pi*f))).^m));
%fircoefs is FIR filter coefficients.  Insert 7 zeros between each:
exp_fircoefs=zeros(1,length(fircoefs)*8-1);
for i=1:length(fircoefs);
    exp_fircoefs(i*8)=fircoefs(i);
end;
[Hfir,w]=freqz(exp_fircoefs,sum(fircoefs),length(cic_resp),'whole');
HFIR=transpose(20*log10(abs(Hfir)));

composite_resp=cic_resp+HFIR;

I think it is easier to work in the frequency domain in this type of problem, especially since you
are after the frequency response as your final answer.



David Joseph Bonnici wrote:

> Thanks for not losing your patience. I am going to make a resume of > what my intentions are and list my two questions in order of priority. > I receive your postings with a latency of 5 hours so if someone has > already made postings giving me solutions please discard these > message. > > Question 1: I need to compute the aggregate frequency response of > three cascaded multirate filters. > I was trying the following but I was in mistake: > What I was trying to do was to get the impulse response of all the > multirate filters and the convolve them. (as suggested on the DDC > datasheet of Xilinx) The problem was that I was not considering the > filters as multirate and so they were being referred to just one > sampling frequency (Fs) and this is not true since Fs changes. > After I have recognized that I had a problem I tried a different > solution: > I've inputted an impulse to just the CIC. After getting the response > (which is also decimated) I've inputted the output waveform into the > PFIR and so on until I get the output from the CFIR. The result I > have obtained does not match with published results of the GSM DDC > from Xilinx. I have tested the all on Matlab R13 which does not have > support of the multirate functions. I have used the following code. > > %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% > clc; %GSM Spectral Mask test so as to check the functionality of the > script with that published in Xilinx DDC > clear; > input = zeros(1000,1); input(1) = 1; % An impulse > > %%%%%%%%%%%%%%%%CIC%%%%%%%%%%%%%%%%%%%%%%%%%%%% > q = quantizer([8 0],'fixed'); % Defines the properties of input. > output_cic = cicdecimate(1, 4, 48, input, q); > > figure; > plot(output_cic); title('Time plot for output due to impluse response > of the CIC'); > freq_y = fft(output_cic,4096); > figure; > l=length(freq_y); > fr=0:1/(l-1)*100e6:100e6; > plot(fr,abs(freq_y)); xlabel('Frequency,Hz'); ylabel('Amplitude'); > title('Frequency Response'); > figure; > plot(fr,20*log10(abs(freq_y))); xlabel('Frequency,Hz'); > ylabel('Amplitude in dB'); title('Frequency Response'); > figure; plot(fr, unwrap(angle(freq_y))); xlabel('Frequency,Hz'); > ylabel('Phase'); title('Phase for y'); > %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% > > %%%%%%%%%%%%%%%PFIR%%%%%%%%%%%%%%%%%%%%%%%% > h = [-98,-679,-2016,-3234,-2537,850,6053,12060,18230,23239,25212, > 23239,18230,12060,6053,850,-2537,-3234,-2016,-679,-98]; > > p =1; > q =2; > output_c = upfirdn(output_cic,h,p,q) > figure; > plot(output_c); title('C section'); > freq_y = fft(output_c,4096); > figure; > l=length(freq_y); > fr=0:1/(l-1)*100e6:100e6; > plot(fr,abs(freq_y)); xlabel('Frequency,Hz'); ylabel('Amplitude'); > title('Frequency Response'); > figure; > plot(fr,20*log10(abs(freq_y))); xlabel('Frequency,Hz'); > ylabel('Amplitude in dB'); title('Frequency Response'); > figure; plot(fr, unwrap(angle(freq_y))); xlabel('Frequency,Hz'); > ylabel('Phase'); title('Phase for y'); > %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% > > %%%%%%%%%%%%%CFIR%%%%%%%%%%%%%%%%%%%%%% > h = [1007,-1853,79,807,1633,423,265,175,-527,-1331,-1454,-1087,-721,-149,1008,1985, > 2164,2005,1483,61,-1756,-3134,-3953,-4016,-2714,134,4003,8361,12934,17009,19565,20353, > 19565,17009,12934,8361,4003,134,-2714,-4016,-3953,-3134,-1756,61,1483,2005,2164,1985,1008,-149, > -721,-1087,-1454,-1331,-527,175,265,423,1633,1807,79,-1853,1007]; > > p =1; > q =2; > output_p = upfirdn(output_c,h,p,q) > figure; > plot(output_p); title('P section'); > freq_y = fft(output_p,4096); > figure; > l=length(freq_y); > fr=0:1/(l-1)*100e6:100e6; > plot(fr,abs(freq_y)); xlabel('Frequency,Hz'); ylabel('Amplitude'); > title('Frequency Response'); > figure; > plot(fr,20*log10(abs(freq_y))-293.1); xlabel('Frequency,Hz'); > ylabel('Amplitude in dB'); title('Frequency Response'); > figure; plot(fr, unwrap(angle(freq_y))); xlabel('Frequency,Hz'); > ylabel('Phase'); title('Phase for y'); > %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% > > You might note that I have been using upfirdn so as to emulate the > polyphase filtering. I do not have a better function so as to emulate > the polyphase filtering. I do not have the latest version of Matlab > so I do not have access to better functions. Still it is not supposed > to make difference since the upsampling ratio is 1 and I am supplying > the filter coefficients myself. > > Now is the problem related with the use of upfirdn or is the > underlying taught completely wrong? What is the correct way of doing > it? > > Question Number 2: > In order of importance Question 1 is more important than this one but > if you can highlight a possible solution I would appreciate it. > Andraka has already suggested a solution but I did not understand what > is he aiming to. > I have used the thing of inserting the zeros between the coefficients > for a different purpose. I have used it so as to get the impulse > response of the CIC referenced to the high sampling frequency. > The code is the following > %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% > clc; > clear; > znum = zeros(193,1); > znum(1) = 1; > znum(49) = -4; > znum(97) = 6; > znum(145) = -4; > znum(193) = 1; > zden = zeros(193,1); > zden(1) = 1; > zden(2) = -4; > zden(3) = 6; > zden(4) = -4; > zden(5) = 1; > > input = zeros(200,1); input(1) = 1; % An impulse > output_cic = filter(znum,zden,input); %the CIC is treated as an IIR > > figure; > plot(output_cic); title('Time plot for output due to impluse response > of the CIC'); > freq_y = fft(output_cic,4096); > figure; > l=length(freq_y); > fr=0:1/(l-1)*100e6:100e6; > plot(fr,abs(freq_y)); xlabel('Frequency,Hz'); ylabel('Amplitude'); > title('Frequency Response'); > figure; > plot(fr,20*log10(abs(freq_y))); xlabel('Frequency,Hz'); > ylabel('Amplitude in dB'); title('Frequency Response'); > figure; plot(fr, unwrap(angle(freq_y))); xlabel('Frequency,Hz'); > ylabel('Phase'); title('Phase for y'); > %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% > The code does its job correctly. > > That's all. Hope you can help me. > Many Thanks again > David Joseph Bonnici
-- --Ray Andraka, P.E. President, the Andraka Consulting Group, Inc. 401/884-7930 Fax 401/884-7950 email ray@andraka.com http://www.andraka.com "They that give up essential liberty to obtain a little temporary safety deserve neither liberty nor safety." -Benjamin Franklin, 1759
Thank's Ray for your help. Your papers regarding the DDC has been of
great help as well.