Multirate Filtering for Digital Signal Processing: MATLAB® Applications, Chapter IV: FIR Filters for Sampling Rate Conversion - Exercises
Contents
Exercise 4.1
Study the polyphase decimation and interpolation by modifying program demo_4_1.
a) Design an optimal linear-phase FIR filter to satisfy the Case b tolerance scheme
of Figure 3.5 (b) for the sampling-rate conversion factor M = L =
4.
The peak passband ripple is ap = 0.1 dB, and the
minimal stopband attenuation amounts to as = 60
dB.
Choose the passband edge frequency ωp < π/M.
The sampling frequency at the input of decimator and at the output of interpolator is 12 kHz.
b) Generate the input signal {x[n]} according to your own choice.
c) Compute four polyphase components for the filter designed in point (a).
d) Use the commutative polyphase structure of Figure 4.8 to decimate-by-4 signal {x[n]}.
To perform decimation, modify Step 4 in program demo_4_1.
e) Use the commutative polyphase configuration of
Figure 4.10 to interpolate-by-4 decimated signal obtained in point
(d). To perform interpolation, modify Step 5 in program demo_4_1.
f) Compute the multiplication rate in decimator and in interpolator, respectively.
g) Compute and plot: filter magnitude response, spectrum of the input signal,
spectrum of the decimated signal, and spectrum of the interpolated signal.
h) Verify the commutative configurations used
above:
(i) Evaluate decimation by using the cascade
connection of FIR filter from point (a) and the factor-of-4
downsampler.
(ii) Evaluate interpolation with the cascade of the factor-of-4
up-sampler and FIR filter from point (a).
(iii) Compute and plot spectra of the decimated signal and that of
the interpolated signal and compare with the results obtained in
point (f).
Filter specification: Case b tolerance scheme (from Figure 3.5 of
the book).
Commutative polyphase structure of a decimator (Figure 4.8 of the book).
Commutative polyphase structure of an interpolator (Figure 4.10 of the book).
Used functions:firpmord, firpm, fir2, freqz, downsample, upsample (Signal Processing Toolbox)
clear all, close all % STEP 1-a),g) % Lowpass FIR filter design M = 4; % Decimation factor L = 4; % Interpolation factor ap = 0.1; as = 60; % Specifying passband/stopband ripples Fp = 0.1; Fs = 2/M-Fp; % Specifying passband/stopband boundary frequencies for Case b tolerance scheme dev = [(10^(ap/20)-1)/(10^(ap/20)+1), 10^(-as/20)]; % Passband and stopband ripples F = [Fp Fs]; % Boundary frequencies A = [1 0]; % Desired amplitudes [Nord,Fo,Ao,W] = firpmord(F,A,dev); % Estimating the filter order Nord = Nord+2; % Increasing the Nord by 2 to provide some margin for Case b tolerance scheme h = firpm(Nord,Fo,Ao,W); % Computing the filter coefficients [H,F1] = freqz(h,1,1024,12000); % Computation of filter frequency response figure (1) subplot(4,1,1), plot(F1,20*log10(abs(H))), title('Figure 1: Spectra (commutative structure)'), legend(' Filter magn. resp.','Location','eastoutside') ylabel('Gain, dB'), axis([0,6000,-80,5]) % STEP 2-b),g) Fx = 12000; % Generating the input signal x n = 0:1000; % Time index ff = [0,0.1,0.5,0.6,0.8,1]; a = [0,0,0.05,0.9,0,0]; x = fir2(1000,ff,a); x = 400*x+cos(2*pi*Fp/2*n)+cos(2*pi*3*Fp/4*n)+randn(size(n)); z=x; [X,F1] = freqz(x,1,1024,12000); % Computing the signal spectrum subplot(4,1,2), plot(F1,abs(X)), legend('Spectr. input sign.','Location','eastoutside'), ylabel('|X(F)|'), axis([0,6000,0,500]) Fy=Fx/M; % STEP 3-c) % Polyphase decomposition E = reshape(h,4,length(h)/4); % STEP 4-d),g) % Polyphase downsampling and filtering x0 = x(1:4:length(x)); x1 = [0,x(4:4:length(x))]; x2 = [0,x(3:4:length(x))]; x3 = [0,x(2:4:length(x))]; y0 = filter(E(1,:),1,x0); y1 = filter(E(2,:),1,x1); y2 = filter(E(3,:),1,x2); y3=filter(E(4,:),1,x3); ydec = y0 + y1 + y2 + y3; % Decimated signal [Ydec,F2] = freqz(ydec,1,512,3000); subplot(4,1,3), plot(F2,abs(Ydec)), legend('Spectr. dec. sign.','Location','eastoutside'), ylabel('|Y_d_e_c(F)|'), axis([0,1500,0,150]) % STEP 5-e),g) x = ydec; % Input to the interpolator % Poliphase filtering and up-sampling u0 = L*filter(E(1,:),1,x); u1 = L*filter(E(2,:),1,x); u2 = L*filter(E(3,:),1,x); u3 = L*filter(E(4,:),1,x); U = [u0;u1;u2;u3]; yint = U(:); % Interpolated signal [Yint,F3] = freqz(yint,1,1024,12000); subplot(4,1,4), plot(F3,abs(Yint)); legend(' Spectr. int. sign.','Location','eastoutside'), xlabel('Frequency, Hz') ylabel('|Y_i_n_t(F)|'), axis([0,6000,0,500]) % STEP 6-h) % Computing the multiplication rate in decimator and in interpolator Rm_dec=(Nord+1)*Fx/M; Rm_int=(Nord+1)*Fx/L; disp('Multiplication rates of commutative structures') disp(['Decimator (mult./second): ', num2str(Rm_dec)]); disp(['Interpolator (mult./second): ', num2str(Rm_int)]); % STEP 7-h) xf = filter(h,1,z); xds = downsample(xf,M); yus = upsample(xds,L); yf = filter(L*h,1,yus); figure (2) subplot(4,1,1), plot(F1,20*log10(abs(H))), title('Figure 2: Spectra (cascade structure)'), legend(' Filter magn. resp.','Location','eastoutside'), ylabel('Gain, dB'), axis([0,6000,-80,5]) subplot(4,1,2), plot(F1,abs(X)), legend('Spectr. input sign.','Location','eastoutside'), ylabel('|X(F)|'), axis([0,6000,0,500]) [Xds,Fds] = freqz(xds,1,512,3000); subplot(4,1,3), plot(Fds,abs(Xds)), legend('Spectr. dec. sign.','Location','eastoutside'), ylabel('|X_d_s(F)|'), axis([0,1500,0,150]) [Yf,Fus] = freqz(yf,1,1024,12000); subplot(4,1,4), plot(Fus,abs(Yf)); legend('Spectr. int. sign.','Location','eastoutside'), xlabel('Frequency, Hz') ylabel('|Y_f(F)|'), axis([0,6000,0,500]) % Comparing Figures 1 and 2, we conclude that commutative structure consisting of 4 polyphase subfilters in parallel is successfully verified. % Both figures display the identical results.
Multiplication rates of commutative structures Decimator (mult./second): 60000 Interpolator (mult./second): 60000
![](Chapter_4_exercises_01.png)
![](Chapter_4_exercises_02.png)
Exercise 4.2
Simulate the process of polyphase decimation and interpolation by modifying program demo_4_2.
a) Design an optimal linear-phase FIR filter to satisfy the Case a tolerance scheme
of Figure 3.5 (a) for the sampling-rate conversion factor M
= L = 6.
The peak passband ripple is ap = 0.1 dB, and the
minimal stopband attenuation amounts to as = 60
dB.
Choose the stopband edge frequency ωs < π/M.
The sampling frequency at the input of decimator and at the output of interpolator is 12 kHz.
b) Generate the input signal {x[n]} according to your own choice.
c) Compute six polyphase components for the filter designed in point (a).
d) Simulate decimation-by-6 of the signal {x[n]} by making use of the commutative polyphase structure
of Figure 4.8. To perform decimation, modify Step 4 in program demo_4_2.
e) Simulate interpolation-by-6 of the decimated signal obtained in point (d) by making use of the commutative
polyphase structure of Figure 4.10. To perform interpolation, modify Step 5 in program demo_4_2.
f) Compute the multiplication rate in decimator and in interpolator, respectively.
g) Compute and plot: filter magnitude response, spectrum of the input signal,
spectrum of the decimated signal, and spectrum of the interpolated signal.
h) Verify the commutative configurations used
above:
(i) Evaluate decimation by using the cascade
connection of FIR filter from point (a) and the factor-of-6
downsampler.
(ii) Evaluate interpolation with the cascade of
the factor-of-6 up-sampler and FIR filter from point (a).
(iii) Compute and plot spectra of the decimated signal and that of the interpolated signal and compare with the results obtained
in point (f).
Filter specification: Case a tolerance scheme (from Figure 3.5 of
the book).
Commutative polyphase structure of a decimator (Figure 4.8 of the book).
Commutative polyphase structure of an interpolator (Figure 4.10 of the book).
Used functions: firpmord, firpm, freqz, (Signal Processing Toolbox)
clear all, close all % STEP 1-a),g) % Lowpass FIR filter design M = 6; % Decimation factor L = 6; % Interpolation factor ap = 0.1; as = 60; % Specifying passband/stopband ripples Fp = 0.1; Fs = 1/M; % Specifying passband/stopband boundary frequencies for Case a tolerance scheme dev = [(10^(ap/20)-1)/(10^(ap/20)+1), 10^(-as/20)]; % Passband and stopband ripples F = [Fp Fs]; % Boundary frequencies A = [1 0]; % Desired amplitudes [Nord,Fo,Ao,W] = firpmord(F,A,dev); % Estimating the filter order Nord = Nord+7; % Increasing the Nord by 7 to provide some margin for Case a tolerance scheme h = firpm(Nord,Fo,Ao,W); % Computing the filter coefficients [H,F1] = freqz(h,1,1024,12000); % Computation of filter frequency response figure (1) subplot(4,1,1), plot(F1,20*log10(abs(H))), title('Figure 1: Spectra (commutative structure)'), legend(' Filter magn. resp.','Location','eastoutside') ylabel('Gain, dB'), axis([0,6000,-80,5]) % STEP 2-b),g) Fx = 12000; % Generating the input signal n = 0:1000; % Time index ff = [0,0.1,0.5,0.6,0.8,1]; a = [0,0,0.05,0.9,0,0]; x = fir2(1000,ff,a); x = 300*x+cos(2*pi*Fp/3*n)+cos(2*pi*Fp/4*n)+0.8*cos(pi*Fs*n)+randn(size(n)); z=x; [X,F1] = freqz(x,1,1024,12000); % Computing the signal spectrum subplot(4,1,2), plot(F1,abs(X)), legend('Spectr. input sign.','Location','eastoutside'), ylabel('|X(F)|'), axis([0,6000,0,500]) % STEP 3-c) % Polyphase decomposition E = reshape(h,M,length(h)/M); % STEP 4-d),g) % Initial states in polyphase components vk_i = zeros(size(E)); % Decimation rot = 0; % Initial rotator position ydec = []; for n=1:length(x) xn = x(n); vk = vk_i(rot+1,:); ek = E(rot+1,:); vk = [xn,vk(1:length(vk)-1)]; yn = sum(vk.*ek); % filtering vk_i(rot+1,:)=vk; if n==1, yy=0; end; if rot == 0, ss = 5; yy = yy+yn; ydec = [ydec,yy]; end; if rot == 1, ss = 0; yy = yy+yn; end; if rot == 2, ss = 1; yy = yy+yn; end; if rot == 3, ss = 2; yy = yy+yn; end; if rot == 4, ss = 3; yy = yy+yn; end; if rot == 5, ss = 4; yy = yn; end; rot = ss; end [Ydec,F2] = freqz(ydec,1,512,2000); subplot(4,1,3), plot(F2,abs(Ydec)), legend('Spectr. dec. sign.','Location','eastoutside') ylabel('|Y_d_e_c(F)|'), axis([0,1000,0,150]) % STEP 5-e),g) % Interpolation yi = 0; xk_i = zeros(size(E)); rot = 0; k = 1; for n=1:length(ydec)*L xn=ydec(k); xk = xk_i(rot+1,:); ek = E(rot+1,:); xk = [xn,xk(1:length(xk)-1)]; yn=sum(xk.*ek); % filtering xk_i(rot+1,:) = xk; yn=L*yn; if rot == 0 ss = 1; yi = [yi,yn]; end; if rot == 1 ss = 2; yi = [yi,yn]; end; if rot == 2 ss = 3; yi = [yi,yn]; end; if rot == 3 ss = 4; yi = [yi,yn]; end; if rot == 4 ss = 5; yi = [yi,yn]; end; if rot == 5 ss = 0; k = k+1; yi = [yi,yn]; end; rot = ss; end [Yint,F3] = freqz(yi,1,1024,12000); subplot(4,1,4), plot(F3,abs(Yint)); legend('Spectr. int. sign.','Location','eastoutside') xlabel('Frequency, Hz'), ylabel('|Y_i_n_t(F)|'), axis([0,6000,0,500]) % STEP 7-f) % Computing the multiplication rate in decimator and in interpolator Rm_dec = (Nord+1)*Fx/M; Rm_int = (Nord+1)*Fx/L; disp('Multiplication rates of commutative structures') disp(['Decimator (mult./second): ', num2str(Rm_dec)]); disp(['Interpolator (mult./second): ', num2str(Rm_int)]); % STEP 8-h) xf = filter(h,1,z); xds = downsample(xf,M); yus = upsample(xds,L); yf = filter(L*h,1,yus); figure (2) subplot(4,1,1), plot(F1,20*log10(abs(H))), title('Figure 2: Spectra (cascade structure)'),legend(' Filter magn. resp.','Location','eastoutside') ylabel('Gain, dB'), axis([0,6000,-80,5]); subplot(4,1,2), plot(F1,abs(X)), legend('Spectr. input sign.','Location','eastoutside'), ylabel('|X(F)|'), axis([0,6000,0,500]) [Xds,Fds] = freqz(xds,1,512,3000); subplot(4,1,3), plot(Fds,abs(Xds)), legend('Spectr. dec. sign.','Location','eastoutside'), ylabel('|X_d_s(F)|'), axis([0,1500,0,150]) [Xus,Fus] = freqz(yf,1,1024,12000); subplot(4,1,4), plot(Fus,abs(Xus)); legend(' Spectr. int. sign.','Location','eastoutside'), xlabel('Frequency, Hz'), ylabel('|Y_f(F)|'), axis([0,6000,0,500]) % Comments % Exercise 4.2 is developed to simulate a sample-by-sample operation in the time-domain for a factor-of-6 decimator and for a factor-of-6 interpolator % based on the efficient commutative polyphase structures. % Comparing the figures 1 and 2 we conclude that commutative structure consisted of 6 polyphase subfilters in parallel is successfully verified. % Both figures display the identical results.
Multiplication rates of commutative structures Decimator (mult./second): 180000 Interpolator (mult./second): 180000
![](Chapter_4_exercises_03.png)
![](Chapter_4_exercises_04.png)
Exercise 4.3
Study the application of memory-saving structures of polyphase decimators
and interpolators by modifying program demo_4_3.
a) Design an optimal linear-phase FIR filter to satisfy the Case a tolerance scheme of Figure 3.5 (a),
for the sampling-rate conversion factor M = L = 4.
The peak passband ripple is ap = 0.1 dB, and the
minimal stopband attenuation is as = 50 dB.
Choose the stopband edge frequency ωs <
π/M.
The sampling frequency at the input of decimator and at the output
of interpolator is 12 kHz.
b) Generate the input signal {x[n]} according to
your own choice.
c) Compute four polyphase components for the
filter designed in point (a).
d) On the basis of the memory-saving structure of Figure 4.13 (from the book) develop the corresponding structure for M = 4.
Decimate-by-4 signal {x[n]} by modifying Step 4 in the program demo_4_3.
e) On the basis of the memory-saving structure of
Figure 4.15 (from the book) develop the corresponding structure for
L = 4. Interpolate-by-4 decimated signal by modifying Step 5 in the program demo_4_3.
f) Compute the multiplication rate in decimator and in interpolator, respectively.
g) Compute the number of delays in decimator and
in interpolator, respectively.
h) Compute and plot: filter magnitude response, spectrum of the input signal,
spectrum of the decimated signal, and spectrum of the interpolated signal.
Filter specification: Case a tolerance scheme (from Figure 3.5 of
the book).
Used functions: firpm, fir2, freqz (Signal Processing Toolbox)
clear all, close all % STEP 1: % Lowpass FIR filter design M = 4; % Decimation factor L = 4; % Interpolation factor Fx = 12000; Fp = 0.15; Fs = 1/M; ap = 0.1; as = 50; % Filter design parameters h = firpm(55,[0,Fp,Fs,1],[1,1,0,0]); % FIR filter design [H,F1] = freqz(h,1,1024,12000); % Computing filter frequency response figure (1) subplot(4,1,1), plot(F1,20*log10(abs(H))), title('Figure 1') legend(' Filter magn. resp.','Location','eastoutside'), ylabel('Gain, dB') axis([0,6000,-60,5]) % STEP 2-b),g) % Generating the input signal {x[n]} n = 0:999; % Time index ff = [0,0.1,0.5,0.6,0.8,1]; a = [0,0,0.05,0.9,0,0]; x = fir2(999,ff,a); x = 300*x+cos(2*pi*Fp/3*n)+cos(2*pi*Fp/4*n)+0.6*cos(pi*Fs*n)+randn(size(n)); z=x; % Input signal [X,F1] = freqz(x,1,1024,12000); % Computing the signal spectrum subplot(4,1,2), plot(F1,abs(X)), legend('Spectr. input sign.','Location','eastoutside'), ylabel('|X(F)|'), axis([0,6000,0,500]) % STEP 3: % Polyphase decomposition E = reshape(h,M,length(h)/M); E = fliplr(E); % STEP 4: % Setting the initial states xk = zeros(size(1:length(h)/M)); % Decimation rot = 0; % Initial rotator position ydec = []; for n=1:length(x) xn = x(n); xk = xn*E((rot+1),:)+xk; if rot == 0 ss = 3; yn = xk(length(xk)); ydec = [yn,ydec]; xk = [0,xk(1:length(xk)-1)]; end; if rot == 1 ss = 0; end; if rot == 2 ss = 1; end; if rot == 3 ss = 2; end; rot = ss; end; [Ydec,F2] = freqz(ydec,1,512,3000); subplot(4,1,3), plot(F2,abs(Ydec)), legend('Spectr. dec. sign.','Location','eastoutside'), ylabel('|Y_{dec}(F)|') axis([0,1500,0,150]) % STEP 5: % Interpolation E = L*fliplr(E); yi = []; xk = zeros(size(1:length(h)/L)); rot = 0; % Initial rotator position k = 1; for n=1:length(ydec)*L xn = ydec(k); ek = E(rot+1,:); if k == 1, xk = [xn,xk(1:length(xk)-1)]; end; yn = sum(xk.*ek); yi = [yi,yn]; if rot == 0 ss = 1; end; if rot == 1 ss = 2; end; if rot == 2 ss = 3; end; if rot == 3 ss = 0; k = k+1; xk = [xn,xk(1:length(xk)-1)]; end; rot = ss; end [Yint,F3] = freqz(yi,1,1024,12000); subplot(4,1,4) plot(F3,abs(Yint)); legend(' Spectr. int. sign.','Location','eastoutside') xlabel('Frequency, Hz'),ylabel('|Y_{int}(F)|') axis([0,6000,0,500]) % STEP 6: % Computation of the multiplication rate in decimator and in interpolator Nf=length(h); % FIR filter length Rm_dec=Nf*Fx/M; % Multiplication rate in decimator Rm_int=Nf*Fx/L; % Multiplication rate in interpolator disp('Multiplication rates of decimator and interpolator') disp(['Decimator (mult./second): ', num2str(Rm_dec)]); disp(['Interpolator (mult./second): ', num2str(Rm_int)]); % STEP 7: % Computation of the number of delays in decimator and in interpolator Nd_dec=Nf/M; % Number of delays in decimator Nd_int=Nf/L; % Number of delays in interpolator disp('Number of delays in decimator and in interpolator') disp(['Decimator: ', num2str(Nd_dec)]); disp(['Interpolator: ', num2str(Nd_int)]); % Comments % The number of delays in the memory-saving structure is determined according to the number of delays in the polyphase component of the highest order. % Filter H(F) of the length N = 56 meets the design requirements specified in Exercise_4_3 and the factor-of-4 decimator and interpolator can be % implemented with 4 polyphase subfilters where the length of each subfilter is 14, and they operate at the intermediate sampling rate of 3 kHz.
Multiplication rates of decimator and interpolator Decimator (mult./second): 168000 Interpolator (mult./second): 168000 Number of delays in decimator and in interpolator Decimator: 14 Interpolator: 14
![](Chapter_4_exercises_05.png)
Exercise 4.4
Investigate the efficiency of two-stage decimator by modifying program .demo_4_4
in such a manner that decimation-by-3 is followed by decimation-by-5.
Generate the corresponding plots to those given in Figures 4.19 and 4.20 (of the book).
Used functions: firpm, freqz, upsample (Signal Processing Toolbox)
close all, clear all h1 = firpm(12,[0,0.5/15,9/15,12/15,1,1],[1,1,0,0,0,0],[1,5,5]); % Designing filter H1(z) [H1,F1] = freqz(h1,1,1024,30); % Frequency response of H1(z) h2 = firpm(59,[0,3*0.5/15,3*1/15,1],[1,1,0,0],[1,5]); % Designing filter H2(z) [H2,F2] = freqz(h2,1,1024,10); % Frequency response of H2(z) h2u = upsample(h2,3); % Periodic filter H_2(z^3) [H2u,F1] = freqz(h2u,1,1024,30); % Spectrum of the periodic filter H_2(z^3) H = H1.*H2u; % The overall frequency response figure (1) subplot(4,1,1), plot(F2,20*log10(abs(H2))), ylabel('|H_2(F)|'), axis([0,5,-80,2]),grid title('Figure 1: Filters: magnitude responses') subplot(4,1,2), plot(F1,20*log10(abs(H2u))), ylabel('|H_2_u(F)|'),axis([0,15,-80,2]),grid subplot(4,1,3), plot(F1,20*log10(abs(H1))), ylabel('|H_1(F)|'), axis([0,15,-80,2]),grid subplot(4,1,4), plot(F1,20*log10(abs(H))), ylabel('|H(F)|'), axis([0,15,-80,2]),grid xlabel('Frequency, kHz') % Input signal ff = [0,1/15,1.5/15,0.6,0.8,1]; a = [0,0,0.8,0.7,0.9,0]; x = fir2(990,ff,a); x = 5*x+0.01*sin(2*pi*(1:991)*300/30000) + 0.01*sin(2*pi*(1:991)*400/30000); [X,F1] = freqz(x,1,2048,30); % Spectrum of the input signal figure (2) subplot(3,1,1), plot(F1,abs(X)), ylabel('|X(F)|'), legend('Original spectrum') title('Figure 2: Signal spectra') % Polyphase decomposition, M1 = 3 h1 = [h1,0,0]; E1 = reshape(h1,3,length(h1)/3); % Polyphase components % Polyphase downsampling and filtering, M1 = 3 v0 = x(1:3:length(x)); v1 = [0,x(3:3:length(x))]; v2 = [0,x(2:3:length(x))]; y0 = filter(E1(1,:),1,v0); y1 = filter(E1(2,:),1,v1); y2 = filter(E1(3,:),1,v2); yd1 = y0 + y1 + y2; % Decimated signal [Yd1,F2]=freqz(yd1,1,2048,10); subplot(3,1,2), plot(F2,abs(Yd1)), legend('1st stage spectrum'), ylabel('|Y_d_1(F)|') % Polyphase decomposition, M2 = 5 E2 = reshape(h2,5,length(h2)/5); % Polyphase downsampling and filtering, M2 = 5 x = yd1; v0 = x(1:5:length(x)); v1 = [0,x(5:5:length(x))]; v2 = [0,x(4:5:length(x))]; v3 = [0,x(3:5:length(x))]; v4 = [0,x(2:5:length(x))]; y0 = filter(E2(1,:),1,v0); y1 = filter(E2(2,:),1,v1); y2 = filter(E2(3,:),1,v2); y3 = filter(E2(4,:),1,v3); y4 = filter(E2(5,:),1,v4); y = y0 + y1 + y2 + y3 + y4; % decimated signal [Y,F3] = freqz(y,1,2048,2); subplot(3,1,3), plot(F3,abs(Y)), legend('2nd stage spectrum'), ylabel('|Y(F)|') xlabel('Frequency, kHz')
![](Chapter_4_exercises_06.png)
![](Chapter_4_exercises_07.png)
Exercise 4.5
Modify program demo_4_4 to study the two-stage
implementation of factor-of-15 interpolator.
a) Implement the interpolator with
L1 = 5 in the first stage, and
L2 = 3 in the second stage.
b) Implement the interpolator with
L1 = 3 in the first stage, and
L2 = 5 in the second stage.
c) Compare the computational efficiency of two solutions from (a) and (b).
Used functions: firpm, freqz, upsample (Signal Processing Toolbox)
close all, clear all % STEP a): Design of two-stage factor-of-15 interpolator with L1=5 and L2=3. L1=5; L2=3; Fp=0.6; ap1=0.01/2; as=60; Fx=2000; % Design parameters for H1(z) Nord_h1=73; % Selecting the filter order for H1(z) h1 = firpm(Nord_h1,[0,0.6/5,1/5,1],[1,1,0,0],[1,3]); % Filter H1(z) [H1,F1] = freqz(h1,1,1024,10); % Frequency response of H1(z) Nord_h2=9; % Selecting the filter order for H2(z) h2 = firpm(Nord_h2,[0,0.6/15,9.4/15,11.6/15,1,1],[1,1,0,0,0,0],[1,3,3]); % Filter H2(z) [H2,F2] = freqz(h2,1,1024,30); % Frequency response of H2(z) h1u = upsample(h1,3); % Periodic filter H_2(z^5) [H1u,F3] = freqz(h1u,1,1024,30); % Spectrum of the periodic filter H_2(z^3) H = H1u.*H2; % The overall frequency response figure (1) subplot(4,1,1), plot(F1,20*log10(abs(H1))), ylabel('|H_1(F)|'), axis([0,5,-80,2]),grid title('Figure 1: Two-stage interpolator L=5x3, equivalent frequency responses, Gain [dB]'); subplot(4,1,2), plot(F3,20*log10(abs(H1u))), ylabel('|H_1_u(F)|'),axis([0,15,-80,2]),grid subplot(4,1,3), plot(F2,20*log10(abs(H2))), ylabel('|H_2(F)|'), axis([0,15,-80,2]),grid subplot(4,1,4), plot(F2,20*log10(abs(H))), ylabel('|H(F)|'), axis([0,15,-80,2]),grid xlabel('Frequency, kHz') % Input signal ff = [0,0.5,0.6,1]; a = [0,0.9,0,0]; x = fir2(100,ff,a); [X,F1] = freqz(x,1,2048,2); % Spectrum of the input signal figure (2) subplot(3,1,1), plot(F1,abs(X)), ylabel('|X(F)|') title('Figure 2: Two stage interpolation, L=5x3: a) Spectrum of the input signal x'); % Polyphase decomposition for the first stage with L1=5 h1=[h1,0]; E1 = reshape(h1,5,length(h1)/5); % Input to the interpolator L1=5 % Poliphase filtering and up-sampling u0 = L1*filter(E1(1,:),1,x); u1 = L1*filter(E1(2,:),1,x); u2 = L1*filter(E1(3,:),1,x); u3 = L1*filter(E1(4,:),1,x); u4 = L1*filter(E1(5,:),1,x); U = [u0;u1;u2;u3;u4]; yint = U(:); % Interpolated signal [Yint,F3] = freqz(yint,1,1024,10); subplot(3,1,2), plot(F3,abs(Yint)); ylabel('|Y_i_n_t(F)|'), axis([0,5,0,5]) title('b) Interpolation-by-5, spectrum of the input signal x after the first stage of interpolation') % Polyphase decomposition for the second stage with L2=3 z=yint'; h2=[h2,0,0]; E2=reshape(h2,3,length(h2)/3); % Input to the interpolator L2=3 % Poliphase filtering and up-sampling v0 = L2*filter(E2(1,:),1,z); v1 = L2*filter(E2(2,:),1,z); v2 = L2*filter(E2(3,:),1,z); V = [v0;v1;v2]; y = V(:); % Interpolated signal [Y,F3] = freqz(y,1,2048,30); subplot(3,1,3), plot(F3,abs(Y)); ylabel('|Y(F)|'), axis([0,15,0,15]) title('c) Interpolation-by-3, spectrum of the input signal x after the second stage of interpolation'), xlabel('Frequency, Hz') % STEP b): Design of two-stage factor-of-15 interpolator with L1=3 and L2=5. L1=3; L2=5; Nord_h3=44; h3 = firpm(Nord_h3,[0,0.6/3,1/3,1],[1,1,0,0],[1,3]); % Filter H3(z) [H3,F1] = freqz(h3,1,1024,10); % Frequency response of H3(z) Nord_h4=16; h4 = firpm(Nord_h4,[0,0.6/15,5.4/15,6.6/15,11.4/15,12.6/15,1,1],[1,1,0,0,0,0,0,0],[1,3,3,3]); % Filter H4(z) [H4,F2] = freqz(h4,1,1024,30); % Frequency response of H4(z) h3u = upsample(h3,5); % Periodic filter H_3(z^5) [H3u,F3] = freqz(h3u,1,1024,30); % Spectrum of the periodic filter H_3(z^5) H = H3u.*H4; % The overall frequency response figure (3) subplot(4,1,1), plot(F1,20*log10(abs(H3))), ylabel('|H_3(F)|'), axis([0,5,-80,2]),grid title('Figure 3: Two-stage interpolator L=3x5, equivalent frequency responses, Gain [dB]'); subplot(4,1,2), plot(F3,20*log10(abs(H3u))), ylabel('|H_3_u(F)|'),axis([0,15,-80,2]),grid subplot(4,1,3), plot(F2,20*log10(abs(H4))), ylabel('|H_4(F)|'), axis([0,15,-80,2]),grid subplot(4,1,4), plot(F2,20*log10(abs(H))), ylabel('|H(F)|'), axis([0,15,-80,2]),grid xlabel('Frequency, kHz') % Input signal: x from Step a). [X,F1] = freqz(x,1,2048,2); % Spectrum of the input signal figure (4) subplot(3,1,1), plot(F1,abs(X)), ylabel('|X(F)|') title('Figure 4: Two stage interpolation, L=3x5: a) Spectrum of the input signal x'); % Polyphase decomposition for the first stage with L1=3 E1 = reshape(h3,3,length(h3)/3); % Input to the interpolator L1=3 % Poliphase filtering and up-sampling u0 = L1*filter(E1(1,:),1,x); u1 = L1*filter(E1(2,:),1,x); u2 = L1*filter(E1(3,:),1,x); U = [u0;u1;u2]; yint = U(:); % interpolated signal [Yint,F3] = freqz(yint,1,1024,6); subplot(3,1,2), plot(F3,abs(Yint)); ylabel('|Y_i_n_t(F)|'), axis([0,3,0,3]) title('b) Interpolation-by-3, spectrum of the input signal x after the first stage of interpolation') % Polyphase decomposition for the second stage with L2=3 z=yint'; h4=[h4,0,0,0]; E2=reshape(h4,5,length(h4)/5); % Input to the interpolator L2=5 % Poliphase filtering and up-sampling v0 = L2*filter(E2(1,:),1,z); v1 = L2*filter(E2(2,:),1,z); v2 = L2*filter(E2(3,:),1,z); v3 = L2*filter(E2(4,:),1,z); v4 = L2*filter(E2(5,:),1,z); V = [v0;v1;v2;v3;v4]; y = V(:); % Interpolated signal [Y,F3] = freqz(y,1,2048,30); subplot(3,1,3), plot(F3,abs(Y)); ylabel('|Y(F)|'), axis([0,15,0,15]) title('c) Interpolation-by-5, spectrum of the input signal x after the second stage of interpolation'), xlabel('Frequency, kHz') % STEP c): Comparation of the computational efficiency of two solutions denoted above as (a) and (b). Nord_h1=73; Nord_h3=44; Nord_h2=9; Nord_h4=16; Fy=15*Fx; RMint_53=(Nord_h1+1)*Fy/15+(Nord_h2+1)*Fy/3; RMint_35=(Nord_h3+1)*Fy/15+(Nord_h4+1)*Fy/5; disp('Multiplication rates of two-stage interpolators:') disp('a) L1 = 5, L2 = 3; b) L1 = 3, L2 = 5') disp(['a) (mult./second): ', num2str(RMint_53)]); disp(['b) (mult./second): ', num2str(RMint_35)]); disp('Comment:') disp('Comparing bottom subfigures of figures 2 and 4 we can conclude that there is no difference in the spectrum shapes of the resulting signals.') disp('Accordingly, both implementations can be considered correct.') disp('Comparing multiplication rates we find that solution b) exhibits better computational complexity.')
Multiplication rates of two-stage interpolators: a) L1 = 5, L2 = 3; b) L1 = 3, L2 = 5 a) (mult./second): 248000 b) (mult./second): 192000 Comment: Comparing bottom subfigures of figures 2 and 4 we can conclude that there is no difference in the spectrum shapes of the resulting signals. Accordingly, both implementations can be considered correct. Comparing multiplication rates we find that solution b) exhibits better computational complexity.
![](Chapter_4_exercises_08.png)
![](Chapter_4_exercises_09.png)
![](Chapter_4_exercises_10.png)
![](Chapter_4_exercises_11.png)
Exercise 4.6
Design and implement the factor-of-6
decimator, which converts the input frequency of 24 kHz to the frequency
of 4 kHz.
a) Design and implement the single-stage
decimator.
b) Design and implement two-stage decimator with
M1 = 3, and M2 = 2.
c) Design and implement two-stage decimator with
M1 = 2, and M2 = 3.
d) Plot and comment on the results.
e) Compute and compare the multiplication rates of the solutions from (a), (b) and (c).
Used functions: firpmord, firpm, fir2, freqz, upsample (Signal Processing Toolbox)
close all, clear all % STEP a) Design and implementation of the factor-of-6 single-stage % decimator. % Lowpass FIR filter design, Case a tolerance scheme. M = 6; % Decimation factor ap = 0.01; as = 60; Fp = 0.05; Fs = 1/M; % FIR filter design parameters Fx = 24000; % Sampling frequency (input) dev = [(10^(ap/20)-1)/(10^(ap/20)+1), 10^(-as/20)]; % Passband and stopband ripples F = [Fp Fs]; % Cutoff frequencies A = [1 0]; % Desired amplitudes [Nord,Fo,Ao,W] = firpmord(F,A,dev); % Estimating the filter order h = firpm(Nord,Fo,Ao,W); % Computing the filter coefficients [H,F1] = freqz(h,1,1024,24); % Computation of filter frequency response figure (1) subplot(3,1,1), plot(F1,20*log10(abs(H))), legend('Filter magnitude response'), ylabel('Gain, dB') title('Figure 1: Single-stage decimator, M=6') axis([0,12,-80,5]),grid Ndec1=Nord; % Generating the input signal n = 0:1008; % Time index ff = [0,0.2,0.5,0.6,0.8,1]; a = [0,0.3,0.5,0.9,0.4,0]; x = fir2(1008,ff,a); x = 300*x+cos(2*pi*Fp/1.9*n)+cos(2*pi*Fp/4*n)+randn(size(n)); z=x; [X,F1] = freqz(x,1,1024,24); % Computing the signal spectrum subplot(3,1,2), plot(F1,abs(X)), legend('Spectrum: input signal'), ylabel('|X(F)|'), axis([0,12,0,500]) Fy=Fx/M; % Sampling frequency (output) % Polyphase decomposition E = reshape(h,6,length(h)/6); % Polyphase downsampling and filtering x0 = x(1:6:length(x)); x1 = [0,x(6:6:length(x))]; x2 = [0,x(5:6:length(x))]; x3 = [0,x(4:6:length(x))]; x4 = [0,x(3:6:length(x))]; x5 = [0,x(2:6:length(x))]; y0 = filter(E(1,:),1,x0); y1 = filter(E(2,:),1,x1); y2 = filter(E(3,:),1,x2); y3=filter(E(4,:),1,x3); y4=filter(E(5,:),1,x4); y5=filter(E(6,:),1,x5); ydec = y0+y1+y2+y3+y4+y5; % Decimated signal [Ydec,F2] = freqz(ydec,1,512,4); subplot(3,1,3), plot(F2,abs(Ydec)), legend('Spectrum: decimated signal'), ylabel('|Y_d_e_c(F)|'), axis([0,2,0,120]) xlabel('Frequency, kHz') % STEP b) Design and implementation of the two-stage decimator M = M1*M2, with M1 = 3, and M2 = 2 M1 = 3; M2 = 2; % Decimation factors ap1 = 0.01/2; as1 = 60; Fp = 0.05; Fs = 1/M1; Fp2 = 3*Fp; Nord_h2 = 20; h2 = firpm(Nord_h2,[0,Fp*3,1/M2,1],[1,1,0,0],[2*W(1),1]); % Second stage, designing filter H2(z) [H2,F2] = freqz(h2,1,1024,8); % Frequency response of H2(z) h2u = upsample(h2,3); % Periodic filter H_2(z^3) [H2u,F1] = freqz(h2u,1,1024,24); % Spectrum of the periodic filter H_2(z^3) Nord_h1 = 11; h1 = firpm(Nord_h1,[0,Fp,6.60/12,9.30/12],[1,1,0,0],[2*W(1),1]); % First stage, designing filter H1(z) [H1,F1] = freqz(h1,1,1024,24); % Frequency response of H1(z) H32 = H1.*H2u; % The overall frequency response figure (2) subplot(4,1,1), plot(F2,20*log10(abs(H2))), legend('Second-stage filter H_2(z)'), ylabel('Gain, dB'), axis([0,4,-80,2]),grid title('Figure 2: Single-stage equivalent of the two-stage decimator, M=3x2') subplot(4,1,2), plot(F1,20*log10(abs(H2u))),legend('Periodic filter H_2(z^3)'), ylabel('Gain, dB'),axis([0,12,-80,2]),grid subplot(4,1,3), plot(F1,20*log10(abs(H1))), legend('First-stage filter H_1(z)'), ylabel('Gain, dB'), axis([0,12,-80,2]),grid subplot(4,1,4), plot(F1,20*log10(abs(H32))), legend('Single-stage equivalent'), ylabel('Gain, dB'),xlabel('Frequency, kHz'), grid % Input signal x1 = z; [X1,F1] = freqz(x1,1,2048,24); % Spectrum of the input signal figure (3) subplot(3,1,1), plot(F1,abs(X1)), ylabel('|X1(F)|'), axis([0,12,0,500]) title('Figure 3: Two stage decimation, M=3x2: Spectrum of the input signal') % Polyphase decomposition, M1 = 3 E1 = reshape(h1,3,length(h1)/3); % Polyphase downsampling and filtering, M1 = 3 v0 = x(1:3:length(x)); v1 = [0,x(3:3:length(x))]; v2 = [0,x(2:3:length(x))]; y0 = filter(E1(1,:),1,v0); y1 = filter(E1(2,:),1,v1); y2 = filter(E1(3,:),1,v2); yd1 = y0 + y1 + y2; % Decimated signal [Yd1,F12] = freqz(yd1,1,2048,8); subplot(3,1,2), plot(F12,abs(Yd1)), ylabel('|Y_d_1(F)|'), title('Decimation-by-3, resulting spectrum after the first stage') % Polyphase decomposition, M2 = 2 h2=[h2,0]; E2 = reshape(h2,2,length(h2)/2); % Polyphase downsampling and filtering, M2 = 2 x = yd1; v0 = x(1:2:length(x)); v1 = [0,x(2:2:length(x))]; y0 = filter(E2(1,:),1,v0); y1 = filter(E2(2,:),1,v1); ydec32 = y0 + y1; % Decimated signal [Ydec32,F3] = freqz(ydec32,1,2048,4); subplot(3,1,3), plot(F3,abs(Ydec32)), ylabel('|Ydec32(F)|'), axis([0,2,0,120]) title('Decimation-by-2, resulting spectrum after the second stage') xlabel('Frequency, kHz') % STEP c) Design and implementation of the two-stage decimator M = M1*M2, with M1 = 2, and M2 = 3 M1 = 2; M2 = 3; % Decimation factors Nord_h3=6; h3 = firpm(Nord_h3,[0,Fp,10.5/12,1],[1,1,0,0],[2*W(1),1]); % Filter H1(z) h3=[h3,0]; [H3,F1] = freqz(h3,1,1024,24); % Frequency response of H1(z) Nord_h4=31; h4 = firpm(Nord_h4,[0,2*Fp,1/M2,1],[1,1,0,0],[2*W(1),1]); % Filter H2(z) [H4,F2] = freqz(h4,1,1024,8); % Frequency response of H2(z) h4u = upsample(h4,2); % Periodic filter H_2(z^2) [H4u,F1] = freqz(h4u,1,1024,24); % Spectrum of the periodic filter H_2(z^2) H23 = H3.*H4u; % The overall frequency response figure (4) subplot(4,1,1), plot(F2,20*log10(abs(H4))), legend('Second-stage filter H4(z)'),ylabel('Gain, dB'), axis([0,4,-80,2]),grid title('Figure 4: Single-stage equivalent of the two stage decimator, M=2x3') subplot(4,1,2), plot(F1,20*log10(abs(H4u))), legend('Periodic filter H4(z^2)'), ylabel('Gain, dB'),axis([0,12,-80,2]),grid subplot(4,1,3), plot(F1,20*log10(abs(H3))), legend('First-stage filter H3(z)'), ylabel('Gain, dB'), axis([0,12,-80,2]),grid subplot(4,1,4), plot(F1,20*log10(abs(H23))), legend('Single-stage equivalent'), ylabel('Gain, dB'), axis([0,12,-80,2]),grid xlabel('Frequency, kHz') % Input signal x2=z; [X2,F1] = freqz(x2,1,1024,24); % Computing the signal spectrum figure (5) subplot(3,1,1), plot(F1,abs(X2)), ylabel('|X2(F)|'), axis([0,12,0,500]) title('Figure 5: Two stage decimation, M=2x3: a) Spectrum of the input signal') % Polyphase decomposition, M1 = 2 E1 = reshape(h3,2,length(h3)/2); % Polyphase downsampling and filtering, M1 = 2 v0 = x2(1:2:length(x2)); v1 = [0,x2(2:2:length(x2))]; y0 = filter(E1(1,:),1,v0); y1 = filter(E1(2,:),1,v1); yd1 = y0 + y1; % Decimated signal [Yd1,F22]=freqz(yd1,1,2048,12); subplot(3,1,2), plot(F22,abs(Yd1)), ylabel('|Y_d_1(F)|'), axis([0,6,0,250]) title('Decimation-by-2, resulting spectrum after the first stage') % Polyphase decomposition, M2 = 3 h4=[h4,0]; E2 = reshape(h4,3,length(h4)/3); % Polyphase downsampling and filtering, M2 = 3 x = yd1; v0 = x(1:3:length(x)); v1 = [0,x(3:3:length(x))]; v2 = [0,x(2:3:length(x))]; y0 = filter(E2(1,:),1,v0); y1 = filter(E2(2,:),1,v1); y2 = filter(E2(3,:),1,v2); ydec23 = y0 + y1 + y2; % Decimated signal [Ydec23,F3] = freqz(ydec23,1,2048,4); subplot(3,1,3), plot(F3,abs(Ydec23)), ylabel('|Ydec23(F)|'), axis([0,2,0,120]) title('Decimation-by-3, resulting spectrum after the second stage') xlabel('Frequency, kHz') % STEP e) % Computing the multiplication rate for decimators designed in STEPS a),b) and c) RM_6_dec1 = (Ndec1+1)*Fx/6; RM1_3_dec2 = (Nord_h1+1)*Fx/3; RM2_2_dec2 = 4*Fx/6; RM_32_dec2 = RM1_3_dec2 + RM2_2_dec2; RM1_2_dec3 = (Nord_h3+1)*Fx/2; RM2_3_dec3 = (Nord_h4+1)*Fx/6; RM_23_dec3 = RM1_2_dec3 + RM2_3_dec3; figure (6) plot(F1,20*log10(abs(H)), F1,20*log10(abs(H32)),F1,20*log10(abs(H23))) xlabel('Frequency, kHz'),ylabel('Gain, dB'),axis([0,12,-120,5]) title('Figure 6: Overall magnitude responses of H(z), H_3_2(z), and H_2_3(z)') legend('H(z)','H_3_2(z)','H_2_3(z)') g(3)=axes('Position',[0.3 0.75 0.40 0.1]); plot(F1(1:60),20*log10(abs(H(1:60))),F1(1:60),20*log10(abs(H32(1:60))),F1(1:60),20*log10(abs(H23(1:60)))) axis([0,0.6,-0.01,0.01]), title('Passband zoom') % Comments % Filters in this example have been designed to satisfy Case a tolerance % scheme. % Figure 6 displays the magnitude responses of the single-stage equivalens % for the three solutions considered in the example. % Figures 1, 3, and 5 illustrate the succesful decimations in the cases of % the single-stage decimator M = 6, and in the cases of the two-stage % decimators M1 = 3, M2 = 2, and M1 = 2, M2 = 3. % Filter complexities: % a) Single-stage decimator M=6: Filter order 59 % Multiplication rate: RM_6_dec1 = 240000 multiplications per second. % b) Two-stage decimator M1=3, M2=2: Filter orders 11 and 20 % Multiplication rate: RM_32_dec2 = 112000 multiplications per second. % c) Two-stage decimator M1=2, M2=3: Filter orders 6 and 31 % Multiplication rate RM_23_dec3 = 212000 multiplications per second. % Evidently, the two-stage solution b) with M1 = 3, M2 = 2 exhibist the best % multiplication rate.
![](Chapter_4_exercises_12.png)
![](Chapter_4_exercises_13.png)
![](Chapter_4_exercises_14.png)
![](Chapter_4_exercises_15.png)
![](Chapter_4_exercises_16.png)
![](Chapter_4_exercises_17.png)
Exercise 4.7
Design and implement the factor-of-6
interpolator, which converts the input frequency of 4 kHz to the frequency
of 24 kHz.
a) Design and implement the single-stage
dinterpolator.
b) Design and implement two-stage interpolator with
L1 = 3, and L2 = 2.
c) Design and implement two-stage decimator with
L1 = 2, and L2 = 3.
d) Plot and comment on the results.
e) Compute and compare the multiplication rates of the solutions from (a), (b) and (c).
Used functions: firpmord, firpm, fir2, freqz, upsample (Signal Processing Toolbox)
close all, clear all % STEP a) Design and implementation of factor-of-6 single-stage interpolator % Lowpass FIR filter design, Case a tolerance scheme. L = 6; % Interpolation factor ap = 0.01; as = 60; Fp = 1.6/12; Fs = 1/L; % FIR filter design parameters Fx = 4000; % Sampling frequency (input) Fy = L*Fx; % Sampling frequency (output) dev = [(10^(ap/20)-1)/(10^(ap/20)+1); 10^(-as/20)]; % Passband and stopband ripples F = [Fp Fs]; % Cutoff frequencies A = [1 0]; % Desired amplitudes [Nord,Fo,Ao,W] = firpmord(F,A,dev); % Estimating the filter order % We optimize the design parameters: Nord=Nord-3; W=[1,1]; h = firpm(Nord,Fo,Ao,W); % Filter H1(z) [H,F1] = freqz(h,1,1024,24); % Frequency response of H1(z) figure (1) subplot(3,1,1), plot(F1,20*log10(abs(H))), ylabel('|H_(F)|'), axis([0,12,-80,2]) title('Figure 1: Single stage interpolator, L=6'); L=6; % Generating the input signal ff = [0,0.7,0.8,1]; a = [0,0.9,0,0]; x = fir2(100,ff,a); [X,F1] = freqz(x,1,1024,4); % Spectrum of the input signal figure (1) subplot(3,1,2), plot(F1,abs(X)), ylabel('|X(F)|'); % Polyphase decomposition for the single stage with L=3 E = reshape(h,6,length(h)/6); % Input to the interpolator L=6 % Poliphase filtering and up-sampling u0 = L*filter(E(1,:),1,x); u1 = L*filter(E(2,:),1,x); u2 = L*filter(E(3,:),1,x); u3 = L*filter(E(4,:),1,x); u4 = L*filter(E(5,:),1,x); u5 = L*filter(E(6,:),1,x); U = [u0;u1;u2;u3;u4;u5]; yint = U(:); % interpolated signal [Yint,F2] = freqz(yint,1,1024,24); subplot(3,1,3), plot(F2,abs(Yint)), ylabel('|Y_i_n_t(F)|'), axis([0,12,0,6]) % STEP b) Design of two-stage factor-of-6 interpolator with L1=3 and L2=2. h1 = firpm(103,[0,1.6/6,2/6,1],[1,1,0,0],[1,1]); % Filter H1(z) h1 = [h1,0]; [H1,F1] = freqz(h1,1,1024,12); % Frequency response of H1(z) h2 = firpm(9,[0,1.6/12,10.5/12,1],[1,1,0,0],[5,1]); % Filter H2(F) [H2,F2] = freqz(h2,1,1024,24); % Frequency response of H2(z) h1u = upsample(h1,2); % Periodic filter H_1(z^2) [H1u,F3] = freqz(h1u,1,1024,24); % Spectrum of the periodic filter H_1(z^2) H32 = H1u.*H2; % The overall frequency response figure (2) subplot(4,1,1), plot(F2,20*log10(abs(H2))), ylabel('|H_2(F)|'), axis([0,12,-80,2]),grid title('Figure 2: Single-stage equivalent of the two-stage interpolator, L=3x2'); subplot(4,1,2), plot(F1,20*log10(abs(H1))), ylabel('|H_1(F)|'), axis([0,6,-80,2]),grid subplot(4,1,3), plot(F3,20*log10(abs(H1u))), ylabel('|H_1_u(F)|'),axis([0,12,-80,2]),grid subplot(4,1,4), plot(F2,20*log10(abs(H32))), ylabel('|H(F)|'), axis([0,12,-80,2]),grid xlabel('Frequency, kHz') L1 = 3; L2 = 2; [X,F1] = freqz(x,1,1024,4); % Spectrum of the input signal figure (3) subplot(3,1,1), plot(F1,abs(X)), ylabel('|X(F)|'); title('Figure 3: Two stage interpolation, L=3x2: Spectrum of the input signal x'); % Polyphase decomposition for the first stage with L1=3 h1=[h1]; E1 = reshape(h1,3,length(h1)/3); % Input to the interpolator L1=3 % Poliphase filtering and up-sampling u0 = L1*filter(E1(1,:),1,x); u1 = L1*filter(E1(2,:),1,x); u2 = L1*filter(E1(3,:),1,x); U = [u0;u1;u2]; yint = U(:); % interpolated signal [Yint,F3] = freqz(yint,1,1024,12); subplot(3,1,2), plot(F3,abs(Yint)), ylabel('|Y_i_n_t(F)|'), axis([0,6,0,3]) title('Interpolation-by-3, resulting spectrum after the first stage') % Polyphase decomposition for the second stage with L2=2 z=yint'; E2=reshape(h2,2,length(h2)/2); % Input to the interpolator L2=2 % Poliphase filtering and up-sampling v0 = L2*filter(E2(1,:),1,z); v1 = L2*filter(E2(2,:),1,z); V = [v0;v1]; y = V(:); % Interpolated signal [Y,F3] = freqz(y,1,2048,24); subplot(3,1,3), plot(F3,abs(Y)); xlabel('Frequency, Hz'), ylabel('|Y_(F)|'), axis([0,12,0,6]) title('Interpolation-by-2, resulting spectrum after the second stage') % STEP c) Design of two-stage factor-of-6 interpolator with L1=2 and L2=3. h3 = firpm(68,[0,1.6/4,2/4,1],[1,1,0,0],[1.5,1]); % Filter H3(z) [H3,F4] = freqz(h3,1,1024,8); % Frequency response of H3(z) h4 = firpm(18,[0,1.6/12,6/12,10/12],[1,1,0,0],[1.5,1]); % Filter H4(z) [H4,F5] = freqz(h4,1,1024,24); % Frequency response of H4(z) h3u = upsample(h3,3); % Periodic filter H_3(z^3) [H3u,F6] = freqz(h3u,1,1024,24); % Spectrum of the periodic filter H_3(z^3) H23 = H3u.*H4; % The overall frequency response figure (4) subplot(4,1,1), plot(F5,20*log10(abs(H4))), ylabel('|H_4(F)|'), axis([0,12,-80,2]),grid title('Figure 4: Single-stage equivalent of the two-stage interpolator, L=2x3'); subplot(4,1,2), plot(F4,20*log10(abs(H3))), ylabel('|H_3(F)|'), axis([0,4,-80,2]),grid subplot(4,1,3), plot(F6,20*log10(abs(H3u))), ylabel('|H_3_u(F)|'),axis([0,12,-80,2]),grid subplot(4,1,4), plot(F5,20*log10(abs(H23))), ylabel('|H(F)|'), axis([0,12,-80,2]),xlabel('Frequency, kHz'), grid L1=2; L2=3; [X,F1] = freqz(x,1,1024,4); % Spectrum of the input signal figure (5) subplot(3,1,1), plot(F1,abs(X)), ylabel('|X(F)|'); title('Figure 5: Two-stage interpolation, L=3x2, spectrum of the input signal x'); % Polyphase decomposition for the first stage with L1=3 h3=[h3,0]; E1 = reshape(h3,2,length(h3)/2); % Input to the interpolator L1=3 % Poliphase filtering and up-sampling u0 = L1*filter(E1(1,:),1,x); u1 = L1*filter(E1(2,:),1,x); U = [u0;u1]; yint = U(:); % Interpolated signal [Yint,F3] = freqz(yint,1,1024,12); subplot(3,1,2), plot(F3,abs(Yint)), ylabel('|Y_i_n_t(F)|'), axis([0,6,0,3]) title('Interpolation-by-2, resulting spectrum after the first stage') % Polyphase decomposition for the second stage with L2=3 z=yint'; h4=[h4,0,0]; E2=reshape(h4,3,length(h4)/3); % Input to the interpolator L2=3 % Poliphase filtering and up-sampling v0 = L2*filter(E2(1,:),1,z); v1 = L2*filter(E2(2,:),1,z); v2 = L2*filter(E2(3,:),1,z); V = [v0;v1;v2]; y = V(:); % Interpolated signal [Y,F3] = freqz(y,1,2048,24); subplot(3,1,3), plot(F3,abs(Y)); xlabel('Frequency, Hz'), ylabel('|Y_(F)|'), axis([0,12,0,6]) title('Interpolation-by-3, resulting spectrum after the second stage') % STEP e) Computation of multiplication rates for three designed interpolators in STEPS 1), 2) and 3) Nord_h = 185; Nord_h1 = 92; Nord_h2 = 5; Nord_h3 = 62; Nord_h4 = 9; Fy = 6*Fx; RMint_6 = (Nord_h+1)*Fy/6; RMint_32 = (Nord_h1+1)*Fy/6+(Nord_h2+1)*Fy/2; RMint_23 = (Nord_h3+1)*Fy/6+(Nord_h4+1)*Fy/3; figure (6) F1 = F1*6; plot(F1,20*log10(abs(H)), F1,20*log10(abs(H32)),F1,20*log10(abs(H23))) xlabel('Frequency, kHz'),ylabel('Gain, dB'),axis([0,12,-120,5]) title('Figure 6: Overall magnitude responses of H(z), H_3_2(z), and H_2_3(z)') legend('H(z)','H_3_2(z)','H_2_3(z)') g(3)=axes('Position',[0.35 0.60 0.40 0.1]); plot(F1(1:150),20*log10(abs(H(1:150))),F1(1:150),20*log10(abs(H32(1:150))),F1(1:150),20*log10(abs(H23(1:150)))) axis([0,1.6,-0.01,0.01]), title('Passband zoom') % Comments % Filters in this example have been designed to satisfy Case a tolerance % scheme. % Figure 6 displays the magnitude responses of the single-stage equivalens % for the three solutions considered in the example. % Figures 1, 3, and 5 illustrate the succesful interpolations in the cases of % the single-stage interpolator L=6, and in the cases of the two-stage % interpolators L1=3, L2=2, and L1=2, L2=3. % Filter complexities: % 1. Single-stage interpolator L=6: Filter order 206 % Multiplication rate: RM_6_int1 = 744000 multiplications per second % 2. Two-stage interpolator M1=3, M2=2: Filter orders 103 and 9, respectively % Multiplication rate: RM_32_int2 = 444000 multiplications per second % 3. Two-stage interpolator L1=2, L2=3: Filter orders 68 and 18,respectively % Multiplication rate RM_23_int3 = 332000 multiplications per second % Evidently, the two-stage solution with L1=2, L2=3 exhibist the best % multiplication rate.
![](Chapter_4_exercises_18.png)
![](Chapter_4_exercises_19.png)
![](Chapter_4_exercises_20.png)
![](Chapter_4_exercises_21.png)
![](Chapter_4_exercises_22.png)
![](Chapter_4_exercises_23.png)