Thank's Ray for your help. Your papers regarding the DDC has been of great help as well.
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
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! =-----
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. > > DavidWhen 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. �����������������������������������������������������������������������
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
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