by Ljiljana Milic
Suplemental material for Chapter IX:
9. Multirate Techniques in Filter Design and Implementation
Table of Contents
9. Multirate Techniques in Filter Design and Implementation
9.1 Multistage Narrowband Filters
Multistage Filtering with the Halfband Decimation and interpolation Filters
Example 9.1
Example 9.2
9.2 Estimation of the Conversion Factor
Example 9.3
9.3 Highpass Multirate Filter
Example 9.4
9.4 Structures Based on Complementary Filters and Multirate Techniques
Example 9.5, and Example 9.6
9.1 Multistage Narrowband Filters
The method is convenient for the lowpass, highpass abd bandpass filters having the bandwidth lower than one fourth of the sampling rate.
Multistage Filtering with the Halfband Decimation and interpolation Filters
Multistage filter for M = 2: Decimator, kernel filter, and interpolator.
Example 9.1
This example demonstrates the contribution of subfilters of the three-stage structure to the overall filter magnitude response and also to the aliasing characteristic.
close all, clear all
Filter specifications
F0 = 2000; % Sampling frequency [Hz]
Fp = 230; % Passband edge frequency [Hz]
Fs = 300; % Stopband edge frequency [Hz]
delta = 0.001; % Peak stopband ripple
Filter design
Nord_k = 50; % Setting the kernel filter order
Designing the kernel filter 
hk = firgr(50,[0,2*Fp/(F0/2),2*Fs/(F0/2),1],[1,1,0,0]);
Designing the decimation filter 
hd = firhalfband('minorder',(2*Fs + Fp)/(3*F0/2),delta);
Interpolation filter 
hi = hd; % Interpolation filter
Frequency responses
[Hk,F] = freqz(upsample(hk,2),1,1024,F0);
[Hd,F] = freqz(hd,1,1024,F0);
figure (1)
plot(F,20*log10(abs(Hk)),'k')
hold on
plot(F,20*log10(abs(Hd)),'b')
ylabel('Gain [dB]')
axis([0,1000,-80,5])
hold on
Hi = Hd;
The overall filter H(z)
H = (Hd.*Hk).*Hi;
figure (2)
plot(F,20*log10(abs(H)),'k')
xlabel('Frequency [Hz]'),ylabel('Gain [dB]')
axis([0,1000,-80,5])
hold on
Computing aliasing
n=1:length(hd);
[Hd_al,F]=freqz(hd,1,2048,F0,'whole');
figure (1)
plot(F(1:1024),20*log10(abs(Hd_al(1025:2048))),'r')
title('H_K(z^2)-black, H_D(z),H_I(z)-blue, H_D(-z)-red' )
Hd_a=Hd_al(1025:2048);
Alias=(Hd_a.*Hk).*Hi;
figure (2)
plot(F(1:1024),20*log10(abs(Alias)),'r')
title('Gain responses of multirate filter')
legend('H(z)','Aliasing')
xlabel('Frequency [Hz]')
disp('END OF EXAMPLE 9.1')
Example 9.2
In this example, we use two identical IIR halfband filters for decimation and interpolation filters 
and , and an eliptic mimimal Q-factors (EMQF) filter for the kernel filter .
and , and an eliptic mimimal Q-factors (EMQF) filter for the kernel filter .close all, clear all
Filter specifications
F0 = 2000; % Sampling frequency [Hz]
Fp = 250; % Passband edge frequency [Hz]
Fs = 300; % Stopband edge frequency [Hz]
delta = 0.001; % Peak stopband filter
Designing the kernel filter
Filter order N = 9
[B,A,beta,Ap,As,Fp,F3db]=emqf_igi(9,1/16+1/32+1/256,0.3);
Designing the decimation filter
Filter order N = 5
[bd,ad,z,p,k] = halfbandiir(5,2*F3db*1000/F0);
Frequency responses
[Hk,F] = freqz(upsample(B,2)/2,upsample(A,2),1024,F0);
[Hd,F] = freqz(bd,ad,1024,F0);
figure (1)
plot(F,20*log10(abs(Hk)),'k',F,20*log10(abs(Hd)),'b')
ylabel('Gain [dB]')
axis([0,1000,-80,5])
hold on
Hi = Hd;
The overall filter H(z)
H = (Hd.*Hk).*Hi;
figure (2)
plot(F,20*log10(abs(H)),'k'),ylabel('Gain [dB]')
axis([0,1000,-80,5])
hold on
Computing aliasing
[Hd_al,F]=freqz(bd,ad,2048,F0,'whole');
figure (1)
figure (1)
plot(F(1:1024),20*log10(abs(Hd_al(1025:2048))),'r')
title('H_K(z^2)-black, H_D(z),H_I(z)-blue, H_D(-z)-red' )
xlabel('Frequency [Hz]')
Hd_a=Hd_al(1025:2048);
Alias=(Hd_a.*Hk).*Hi;
figure (2)
plot(F(1:1024),20*log10(abs(Alias)),'r')
title('Gain responses of multirate filter')
xlabel('Frequency [Hz]')
disp('END OF EXAMPLE 9.2')
9.2 Estimation of the Conversion Factor
The optimal conversion factor M is determined as the integer part of the real root of of the equation:
Example 9.3
Specifications for FIR lowpass filter:
Bandpass edge frequency: 
Stopband edge frequency: 
Sampling frequency: 
Ripple tolerance in the passband: 
Minimal stopband attenuation 60 dB (
)
)close all, clear all
F0 = 2000; Fp = 50; Fs = 100;
Estimating the optimal conversion factor
D = [Fs^2-Fp^2,-(Fs+Fp)^2,2*F0*(Fs+Fp),-F0^2];
D_roots = roots(D);
M = fix(D_roots(3));
disp(' EXAMPLE 9.3')
disp('Roots of the polynomial D')
disp(' D_roots = ')
disp(D_roots')
disp('Optimal conversion factor')
disp([sprintf(' M = %15.9f', M)])
Kernel filter design
fk = [0,50/200,100/200,1]; mk=[1,1,0,0]; % Setting the design parameters
w = [1,3.333]; % Passband/stopband weighting coefficients
hk = firpm(24,fk,mk,w); % Computation of filter coefficients
[Hk,fk]=freqz(hk,1,512,400);
Decimation/ interpolation filter design
f0 = [0,50/1000,350/1000,450/1000,750/1000,850/1000]; % Design parameters
m0 = [1,1,0,0,0,0]; % Design parameters
w0 = [w,w(2)]; % Passband/stopband weighting coefficients
hd = firpm(18,f0,m0,w0); % Computation of filter coefficients
[Hd,f]=freqz(hd,1,512,2000); % Decimation/interpolation filter
Interpolated kernel filter
hkm=zeros(size(1:M*length(hk)));
hkm([1:5:length(hkm)])=hk;
[Hkm,f]=freqz(hkm,1,512,2000);
Multirate filter frequency response
H=(Hd.*Hkm).*Hd;
Displaying the results
figure (1)
plot(fk,20*log10(abs(Hk)));
xlabel('Frequency [Hz]'),ylabel('Gain [dB]')
title('Kernel filter H_K(z)')
axis([0,200,-80,5]), grid
figure (2)
plot(f,20*log10(abs(Hd)),'r');
xlabel('Frequency [Hz]'),ylabel('Gain [dB]')
title('Decimation and interpolation filter')
axis([0,1000,-80,5]), grid
figure (3)
plot(f,20*log10(abs(H)));
xlabel('Frequency [Hz]'),ylabel('Gain [dB]')
title(' Multistage filter, H(z)')
axis([0,1000,-90,5]), grid
Interpolation filter
h=hd; % Impulse response of interpolation filter
n=0:length(h)-1;
Fx=2000; % Input sampling frequency
M=5; % Decimation factor
Fy=Fx/M; % Ooutput sampling frequency
for r=1:1001
F(r)=(r-1)*Fx/400;
omega=2*pi*F(r)/Fx;
W=exp(-i*n*omega);
H(r)=sum(h.*W); % FIR filter frequency response
end
Computing and displaying aliasing
Computation of aliasing in the frequency band of the low-rate signal
For k =0 program computes unalised component
For k = 1, 2, 3, 4 program computes the aliased components
k=0;
for l=1:5
for r=1:201
FF(r)=(r-1)*Fy/75;
omega=2*pi*FF(r)/Fy;
W2=exp(-i*n*(omega-2*k*pi)/M);
Hal(r)=sum(h.*W2);
end
k=k+1;
plot(FF,20*log10(abs(Hal)))
hold on
end
xlabel('Frequency [Hz]'),ylabel('Gain [dB]')
title('Aliasing characteristics')
legend('Location','south'), legend('Aliasing')
axis([0,50,-90,5]), grid on
hold off
disp('END OF EXAMPLE 9.3')
9.3 Highpass Multirate Filter
Example 9.4
In this example, we show how the lowpass filter of Example 9.1 is modified into a highpass multirate filter.
Filter specifications
close all, clear all
F0 = 2000;
Fp = 230; % LP filter
Fs = 300; % LP filter
FP_h = 1000-Fp; % HP filter
FS_h = 1000-Fs; % HP filter
delta = 0.001;
Kernel filter design
hk = firgr(50,[0,2*Fp/(F0/2),2*Fs/(F0/2),1],[1,1,0,0]);
Design of highpass halfband decimation filter
hd = firhalfband('minorder',(2*Fs+Fp)/(3*F0/2),delta);
n=0:length(hd)-1;
hd=hd.*(-1).^n;
Frequency responses
[Hk,F] = freqz(upsample(hk,2),1,1024,F0);
[Hd,F] = freqz(hd,1,1024,F0);
figure (1)
plot(F,20*log10(abs(Hk)),F,20*log10(abs(Hd)),'g')
ylabel('Gain [dB]')
axis([0,1000,-80,5])
hold on
Highpass interpolation filter
hi = hd;
Hi = Hd;
Multirate filter frequency response
H = (Hd.*Hk).*Hi;
figure (2)
plot(F,20*log10(abs(H))),ylabel('Gain [dB]')
axis([0,1000,-80,5])
hold on
Computing aliasing and displaying results
n=1:length(hd);
[Hd_al,F]=freqz(hd,1,2048,F0,'whole');
figure (1)
plot(F(1:1024),20*log10(abs(Hd_al(1025:2048))),'r')
legend('location','best')
legend('H_K(z^2)','H_D(z),H_I(z)','H_D(-z)')
xlabel('Frequency (Hz)')
Hd_a=Hd_al(1025:2048);
Alias=(Hd_a.*Hk).*Hi;
figure (2)
plot(F(1:1024),20*log10(abs(Alias)),'r'),grid
title('Multirate filter')
legend('Location','best')
legend('H(z)','aliasing')
xlabel('Frequency [Hz]')
disp('END OF EXAMPLE 9.4')
9.4 Structures Based on Complementary Filters and Multirate Techniques
The following examples show the application of multirate approach to design broadband filters. This is achieved by subtracting output of the narrow-band filter from the delayed version of the input signal.
Example 9.5, and Example 9.6
Example 9.5
In this example, we demonstrate the application of the multirate complementary structure, shown bellow, to construct a wideband lowpass filter. The multirate structure consists of the kernel filter , decimation filter , interpolation filter , and delay element 
.
, decimation filter , interpolation filter , and delay element .
Wideband filter specifications:
Sanpling frwquency:
Passband edge frequency 
Stopband edge frequency 
Minimal stopband attenuation 
The design procedure consists of two steps:
(1) Design the narrowband multirate highpass filter using the procedure of Example 9.4.
(2) Determine the the desired wideband filter as delay-complementary filter to the highpass narrowband filter.
close all, clear all
(1) Designing the multirate highpass filter from Example 9.4
F0 = 2000;
Fp = 230;
Fs = 300;
delta = 0.0003;
Kernel filter design
hk = firgr(62,[0,2*Fp/(F0/2),2*Fs/(F0/2),1],[1,1,0,0]);
Decimation/interpolation filter
hd = firhalfband('minorder',(2*Fs+Fp)/(3*F0/2),delta);
n=0:length(hd)-1;
hd=hd.*(-1).^n;
hi = hd;
Frequency responses
[Hk,F] = freqz(upsample(hk,2),1,1024,F0);
[Hd,F] = freqz(hd,1,1024,F0);
Hi = Hd;
Multirate highpass filter
H = (Hd.*Hk).*Hi;
Computing aliasing
n=1:length(hd);
[Hd_al,F1]=freqz(hd,1,2048,F0,'whole');
figure (1)
Hd_a=Hd_al(1025:2048);
Alias=(Hd_a.*Hk).*Hi;
Displaying highpass filter characteristics
figure (1)
plot(F,20*log10(abs(H))),ylabel('Gain [dB]')
axis([0,1000,-100,5])
hold on
plot(F1(1:1024),20*log10(abs(Alias)),'--b'),grid
title('Multirate highpass filter')
legend('Location','best')
legend('H(z)','aliasing')
xlabel('Frequency [Hz]')
(2) Designing the wideband filter
D=length(hk)+length(hd)-2; % Computing delay D (samples)
Computing the wideband filter frequency response
HH=ones(size(1:length(H))).*exp(-i*2*pi*F*D/F0)'-H';
figure (2)
plot(F,20*log10(abs(HH)),F1(1:1024),20*log10(abs(Alias)),'--b')
ylabel('Gain [dB]'),xlabel('Frequency [Hz]')
title('Multirate wideband filter')
legend('Location','northwest')
legend('Filter','aliasing')
axis([0,1000,-120,5]),grid
g(1)=axes('Position',[0.25 0.50 0.45 0.1]);
plot(F(1:750),(abs(HH(1:750)))), axis([0,730,0.9998,1.0002]),grid
disp('END OF EXAMPLE 9.5')
Example 9.6
This example demonstrates the application of the two-stage complementary structure, shown in the Figure bellow, to construct a sharp lowpass filter with the design parameters:
Sampling frequency: 
Passband edge frequency
Stopband edge frequency
Minimal stopband attenuation 
Solution
For the first step, we exploit the filters from the second stage of Example 9.5. Here, this filter has a role of the equivalent kernel filter. In the second step, we design decimation and interpolation filters 
and 
.
and .disp(' EXAMPLE 9.6')
Computing the impulse response of the equivalent kernel filter:
LP multirate filter of the first stage
hy=conv(downsample(hd,2),hk);
hy=conv(upsample(hy,2),2*hi);
hdel=zeros(size(1:length(hy)));
D=(length(hdel))/2;
hdel(D)=1;
hh=hdel-hy; % Impulse response of the equivalent kernel filter
[HH,F]=freqz(hh,1,1024,F0); % Frequency response
figure (3)
plot(F,20*log10(abs(HH)),F1(1:1024),20*log10(abs(Alias)),'--b')
title('LP multirate filter from the first stage')
ylabel('Gain [dB]'),xlabel('Frequency [Hz]')
grid
axis([0,1000,-120,5])
g(2)=axes('Position',[0.25 0.70 0.45 0.1]);
plot(F(1:750),(abs(HH(1:750))))
axis([0,730,0.9998,1.0002]),grid
Decimation/interpolation filter
hd1 = firhalfband('minorder',(2*Fs+2.2*Fp)/(3*F0/2),delta);
n=0:length(hd1)-1;
hd1_a=hd1.*(-1).^n;
F0=4000; % Sampling frequency of the overall multirate filter
hi1 = hd1;
Computing the impulse response of the overall multirate filter
hy1=conv(downsample(hd1,2),hh);
hy1=conv(upsample(hy1,2),2*hi1); % The output sequence
hdel1=zeros(size(1:length(hy1)));
D1=(length(hdel1))/2;
hdel(D1)=1;
hh1=hdel1-hy1; % Overall impulse response
Frequency responses
[Hd1,F1]=freqz(hd1,1,1024,4000);
[Hd1_a,F1]=freqz(hd1_a,1,1024,4000); % Aliasing
[HHu,F]=freqz(upsample(hh,2),1,1024,4000); % Upsampled eq. kernel filter
[HH1,F]=freqz(hh1,1,1024,F0);
Alias_1=(Hd1_a.*HHu).*Hd1;
Alias_1a=Alias.*Hd1;
Displaying the gain responses of the two-stage wideband multirate filter
figure (4)
plot(F,20*log10(abs(HH1)));
hold on
plot(F,20*log10(abs(Alias_1)),'--b');
plot(F,20*log10(abs(Alias_1a)),'--b')
ylabel('Gain [dB]'),xlabel('Frequency [Hz]')
grid
title('Two-stage multirate LP filter')
legend('Location','best')
legend('LP filter','Aliasing')
axis([0,2000,-120,5])
g(3)=axes('Position',[0.52 0.65 0.37 0.1]);
plot(F(1:450),(abs(HH1(1:450)))), axis([0,730,0.999,1.001]),grid
disp('END OF EXAMPLE 9.6')
disp(' END OF CHAPTER IX')