by Ljiljana Milic
Supplemental material for Chapter III:
3. Filters in Multirate Systems
Table of Contents
3. Filters in Multirate Systems
3.1 Introduction
3.2 FIR Filter Design
Example 3.1
Example 3.2
Example 3.3
Example 3.4
Example 3.5
3.3 IIR Filter Design
Example 3.6
Example 3.7
3.4 Computation of Aliasing Characteristics
Example 3.8
3.5 Sampling Rate Alteration of Bandpass Signals
3.5.1 Bandpass Signal Decimation
Example 3.9
3.5.2 Integer Band Interpolation
Example 3.10
3.1 Introduction
In sampling-rate conversion systems filters are used to suppress aliasing in decimation, and to remove imaging in interpolation. Three toleranc schemes given below show three types of filter specifications for sampling-rate conversion by factor-of-M.
(a) Case a tolerance scheme, (b) Case b tolerance scheme, (c) Case c tolerance scheme. Note that for IIR filterdesign, the passpand magnitude response is limiter to the range 
.
.3.2 FIR Filter Design
Example 3.1
Linear-phase FIR filter designed to satisfy Case a tolerance scheme for the sampling-rate conversion factor M = 5 .
clear all, close all
Filter specifications
ap = 0.1; as = 60; Fp = 0.09; Fs = 1/5; %
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
Estimating the filter order
[Nord,Fo,Ao,W] = firpmord(F,A,dev);
Computing the filter coefficients.
h = firpm(Nord+4,Fo,Ao,W);
Compute and plot frequency response.
[H,f] = freqz(h,1,1024,2);
figure
plot(f,20*log10(abs(H))),grid
xlabel('\omega/\pi')
ylabel('Magnitude, dB')
axis([0,1,-80,5])
k=fix(1024*0.09);
g=axes('Position',[0.40 0.71 0.45 0.1]);
plot(f(1:92),20*log10(abs(H(1:92)))),grid
figure
Phase = unwrap(angle(H));
plot(f,Phase), grid
xlabel('\omega/\pi'),ylabel('Phase, rad')
Plot the impulse response and pole-zero locations.
figure
subplot(2,1,1)
impz(h,1); % Plots the impulse response
subplot(2,1,2)
zplane(h,1) % Pole-zero locations
disp('END OF EXAMPLE 3.1')
Example 3.2
Linear-phase FIR filter designed to satisfy Case b tolerance scheme for the sampling-rate conversion factor M = 5 .
clear all, close all
ap = 0.1; as = 60; Fp = 0.09; Fs = 2/5-Fp; % Filter specifications
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
[N,Fo,Ao,W] = firpmord(F,A,dev); % Estimating the filter order
h = firpm(N+3,Fo,Ao,W);
figure
[H,f] = freqz(h,1,1024,2); % Computes the frequency response
figure
plot(f,20*log10(abs(H))),grid
xlabel('\omega/\pi')
ylabel('Magnitude, dB')
axis([0,1,-80,5])
k=fix(1024*0.09);
g=axes('Position',[0.40 0.71 0.45 0.1]);
plot(f(1:92),20*log10(abs(H(1:92)))),grid
figure
Phase = unwrap(angle(H));
plot(f,Phase), grid
xlabel('\omega/\pi'),ylabel('Phase, rad')
figure
subplot(2,1,1)
impz(h,1) % Plots the impulse response
subplot(2,1,2)
zplane(h,1) % Pole-zero locations
disp('END OF EXAMPLE 3.2')
Example 3.3
Linear-phae FIR filter designed to satisfy Case c tolerance scheme for the sampling-rate conversion factor M = 5 .
clear all, close all
ap=0.1; as=60; M = 5;
Fp=0.09; ;
dev = [(10^(ap/20)-1)/(10^(ap/20)+1),10^(-as/20),10^(-as/20),10^(-as/20)];
F = [Fp,2/M-Fp,2/M+Fp,4/M-Fp,4/M+Fp,1]; % Cutoff frequencies
A = [1,0,0,0]; % Desired amplitudes
[Nord,Fo,Ao,W] = firpmord(F,A,dev);
Nord = 26;
Fop = [0, Fp]; % Passband boundary frequencies
Fos1 = [2/M-Fp,2/M+Fp]; % Stopband 1 boundary frequencies
Fos2 = [4/M-Fp,4/M+Fp]; % Stopband 2 boundary frequencies
Fo= [Fop, Fos1, Fos2]; % Vector of normalized frequency points
Ao = [1,1,0,0,0,0]; W =[1.0000,5.7564,5.7564];
h = firpm(Nord,Fo,Ao,W);
figure
[H,f] = freqz(h,1,1024,2); % Computes the frequency response
figure
plot(f,20*log10(abs(H))),grid
xlabel('\omega/\pi')
ylabel('Magnitude, dB')
axis([0,1,-80,5])
k=fix(1024*0.09);
g=axes('Position',[0.40 0.71 0.45 0.1]);
plot(f(1:92),20*log10(abs(H(1:92)))),grid
figure
Phase = unwrap(angle(H));
plot(f,Phase), grid
xlabel('\omega/\pi'),ylabel('Phase, rad')
figure
subplot(2,1,1)
impz(h,1) % Plots the impulse response
subplot(2,1,2)
zplane(h,1) % Pole-zero plot
disp('END OF EXAMPLE 3.3')
Example 3.4
Minimum-phase FIR filter designed to satisfy Case a tolerance scheme for the sampling-rate conversion factor M = 5 .
clear all, close all
Fp=0.09; ; Fs=1/5;
Nord = 38; Fo = [ 0, Fp, Fs, 1.0000]; Ao = [1, 1, 0, 0]; W = [1.0000, 200];
h = firgr (Nord,Fo,Ao,W,'minphase'); % Filter coefficients
figure
[H,f] = freqz(h,1,1024,2); % Computes the frequency response
plot(f,20*log10(abs(H))),grid
xlabel('\omega/\pi')
ylabel('Magnitude, dB')
axis([0,1,-80,5])
k=fix(1024*0.09);
g=axes('Position',[0.40 0.71 0.45 0.1]);
plot(f(1:92),20*log10(abs(H(1:92)))),grid
figure
Phase = unwrap(angle(H));
plot(f,Phase), grid
xlabel('\omega/\pi'),ylabel('Phase, rad')
figure
subplot(2,1,1)
impz(h,1) % Plots the impulse response
subplot(2,1,2)
zplane(h,1) % Pole-zero plot
disp('END OF EXAMPLE 3.4')
Example 3.5
The interpolated FIR (IFIR) is based on the cascade of two filters where the first filter in cascade is the periodic filter 
, and the second filter 
is is called the image suppression filter. The realization structure of the IFIR filter is depicted bellow, and in the overall transfer function is given by
, and the second filter is is called the image suppression filter. The realization structure of the IFIR filter is depicted bellow, and in the overall transfer function is given by
Design IFIR filter to satisfy Case a tolerance scheme.
Design specifications:
- Sampling-rate conversion factor M = 5.
- Passband edge
. Stopband edge at 
. - Pasband and stopband ripples are given in decibels with
, and 
.
clear all, close all
Filter specifications
L=3; % Choice of the upsampling factor for the periodic filter
ap = 0.1; as = 60;
Fp = 0.09; Fs = 0.2;
dev = [(10^(ap/20)-1)/(10^(ap/20)+1), 10^(-as/20)];
Design of periodic filter and 
and [h,g] = ifir(L,'low',[Fp, Fs],dev);
Computing frequency responses of and
and [H,f] = freqz(h,1,1024,2); [G,f]=freqz(g,1,1024,2);
Computing the overall impulse response of IFIR filter
hIFIR = conv(h,g);
Computing the frequency response of IFIR filter
[HIFIR,f] = freqz(hIFIR,1,1024,2);
Displaying the resulta
figure
subplot(2,1,1)
plot(f,20*log10(abs(H)));
hold on
plot(f,20*log10(abs(G))); grid
axis([0,1,-80,5]), xlabel('\omega/\pi')
ylabel('Magnitude, dB')
legend('Periodic Filter','Image Suppressor Filter');
subplot(2,1,2)
plot(f,20*log10(abs(HIFIR))),grid
axis([0,1,-150,5]), legend('Overall Filter');
xlabel('\omega/\pi'),ylabel('Magnitude, dB')
figure
subplot(2,2,1), stem(0:length(h)-1,h)
ylabel('h[n]'),xlabel('n'), axis([0,51,-0.1,0.4])
subplot(2,2,2), stem(0:length(g)-1,g)
ylabel('g[n]'),xlabel('n'), axis([0,14,-0.1,0.3])
subplot(2,2,3), stem(0:length(hIFIR)-1,hIFIR)
ylabel('h_I_F_I_R[n]'),xlabel('n'), axis([0,length(hIFIR)-1,-0.05,0.15])
subplot(2,2,4), zplane(hIFIR,1)
disp('END OF EXAMPLE 3.5')
3.3 IIR Filter Design
Example 3.6
Elliptic IIR filter designed to satisfy Case a tolerance scheme for the sampling-rate conversion factor M = 5 .
clear all, close all
Design speciffications
ap=0.1; as=60;
Fp=0.09;
N = 5; % Filter order
Elliptic filter design
[B,A]=ellip(5,ap,as,Fp);
Computing the filter frequency response
[H,f] = freqz(B,A,1024,2);
Disply the results
figure
plot(f,20*log10(abs(H))),grid
xlabel('\omega/\pi')
ylabel('Magnitude, dB')
axis([0,1,-80,5])
k=fix(1024*0.09);
g=axes('Position',[0.40 0.71 0.45 0.1]);
plot(f(1:k),20*log10(abs(H(1:92)))),grid
figure
Phase = unwrap(angle(H));
plot(f,Phase), grid
xlabel('\omega/\pi'),ylabel('Phase, rad')
figure
subplot(2,2,1:2)
impz(B,A) % Plots the impulse response
axis([0,60,-0.05,0.15])
subplot(2,2,4)
zplane(B,A) % Pole-zero locations
disp('END OF EXAMPLE 3.6')
Example 3.7
Chebyshev IIR filter designed to satisfy Case a tolerance scheme for the sampling-rate conversion factor M = 5 .
clear all, close all
Design specifications
ap=0.1; as=60;
Fp=0.09;
N = 7; % Filter order
Chabyshev filter design
[B,A]=cheby1(7,ap,Fp);
Computing the filter frequency response
[H,f] = freqz(B,A,1024,2);
Disply the results
figure
plot(f,20*log10(abs(H))),grid
xlabel('\omega/\pi')
ylabel('Magnitude, dB')
axis([0,1,-80,5])
k=fix(1024*0.09);
g=axes('Position',[0.40 0.71 0.45 0.1]);
plot(f(1:k),20*log10(abs(H(1:92)))),grid
figure
Phase = unwrap(angle(H));
plot(f,Phase), grid
xlabel('\omega/\pi'),ylabel('Phase, rad')
figure (2)
subplot(2,2,1:2)
impz(B,A) % Plots the impulse response
axis([0,60,-0.05,0.15])
subplot(2,2,4)
zplane(B,A) % Pole/zero plot
disp('END OF EXAMPLE 3.7')
3.4 Computation of Aliasing Characteristics
In this section, we demonstrate by means of example 3.8 the filter contribution to the suppression of aliasing in decimation.
The role of decimation (antialiasing) filter 
is to attenuate the spectral components of the high-rate signal above the half of the sampling-frequency of the high-rate signal.
is to attenuate the spectral components of the high-rate signal above the half of the sampling-frequency of the high-rate signal.
The spectrum of the lo-rate signal 
is expressible in terms of the spectrum of the high-rate signal 
and the frequency response of the antialiasing filter 
is expressible in terms of the spectrum of the high-rate signal and the frequency response of the antialiasing filter 
where
Here, for k = 0 function 
represents the unalized characteristic, whereas for k = 1,..., M - 1, represents M - 1 aliasing characteristics for a factor-of-M decimator.
represents the unalized characteristic, whereas for k = 1,..., M - 1, represents M - 1 aliasing characteristics for a factor-of-M decimator.To determine the suppression of the aliased spectra, we compute the absolute values of and usually express in decibels,
and usually express in decibels,
For more details see the book chapter III, pages 89 - 93.
Example 3.8
In this example, we compute and and plot aliasing components for FIR antialiasing filters from examples 3.1, 3.2, and 3.3 which were designed to satisfy Case a,Case b, and Case c tolerance schemes, respectively.
Design specifications:
close all, clear all
M = 5; % Decimation factor
ap=0.1; as=60; Fp=0.09; % Specified for cases (a), (b), and (c).
dev = [(10^(ap/20)-1)/(10^(ap/20)+1) 10^(-as/20)];
% Stop-band cutoffs for Case (a), Case (b), and Case (c), respectively:
Fsa=1/M; Fsb = 2/M-Fp; Fsc = [2/M-Fp,2/M+Fp,4/M-Fp,4/M+Fp,1];
Fx = 10000; % Imput sampling frequency
Fy =2000; % Output sampling frequency
Case a design
F = [Fp Fsa]; % Cutoff frequencies
A = [1 0]; % Desired amplitudes
[Nord,Fo,Ao,W] = firpmord(F,A,dev);
ha = firpm(Nord+4,Fo,Ao,W); % Filter coefficients
Computing frequency response for Case a design
n=0:length(ha)-1;
for r=1:1001
F(r)=(r-1)*Fx/2000;
omega=2*pi*F(r)/Fx;
W=exp(-i*n*omega);
Ha(r)=sum(ha.*W);
end
figure
subplot(2,1,1)
plot(F,20*log10(abs(Ha)))
xlabel('Frequency, Hz'),ylabel('Magnitude, dB')
title('Magnitude response for Case a design')
axis([0,5000,-80,5]), grid
Computing aliasing characteristic for Case a design
subplot(2,1,2)
k=0;
for l=1:5
for r=1:201
FF(r)=(r-1)*Fy/400;
omega=2*pi*FF(r)/Fy;
W2=exp(-i*n*(omega-2*k*pi)/M);
Hala(r)=sum(ha.*W2);
end
k=k+1;
plot(FF,20*log10(abs(Hala)))
hold on
end
xlabel('Frequency, Hz'),ylabel('Magnitude, dB')
title('Aliasing for Case a design')
axis([0,1000,-80,5]), grid on
hold off
Case b design
F = [Fp Fsb]; % Cutoff frequencies
A = [1 0]; % Desired amplitudes
[Nord,Fo,Ao,W] = firpmord(F,A,dev);
hb = firpm(Nord+3,Fo,Ao,W); % Filter coefficients
Computing frequency response for Case b design
n=0:length(hb)-1;
for r=1:1001
F(r)=(r-1)*Fx/2000;
omega=2*pi*F(r)/Fx;
W=exp(-i*n*omega);
Hb(r)=sum(hb.*W);
end
figure
subplot(2,1,1)
plot(F,20*log10(abs(Hb)))
xlabel('Frequency, Hz'),ylabel('Magnitude, dB')
title('Magnitude response for Case b design')
axis([0,5000,-80,5]), grid
Computing aliasing characteristic for Case b design
subplot(2,1,2)
k=0;
for l=1:5
for r=1:201
FF(r)=(r-1)*Fy/400;
omega=2*pi*FF(r)/Fy;
W2=exp(-i*n*(omega-2*k*pi)/M);
Halb(r)=sum(hb.*W2);
end
k=k+1;
plot(FF,20*log10(abs(Halb)))
hold on
end
xlabel('Frequency, Hz'),ylabel('Magnitude, dB')
title('Aliasing for Case b design')
axis([0,1000,-80,5]), grid on
hold off
Case c design
F = [Fp Fsc]; % Cutoff frequencies
A = [1,0,0,0]; % Desired amplitudes
dev = [(10^(ap/20)-1)/(10^(ap/20)+1),10^(-as/20),10^(-as/20),10^(-as/20)];
[Nord,Fo,Ao,W] = firpmord(F,A,dev);
hc = firpm(26,Fo,Ao,W); % Filter coefficients
Computing frequency response for Case (c) design
n=0:length(hc)-1;
for r=1:1001
F(r)=(r-1)*Fx/2000;
omega=2*pi*F(r)/Fx;
W=exp(-i*n*omega);
Hc(r)=sum(hc.*W);
end
figure
subplot(2,1,1)
plot(F,20*log10(abs(Hc)))
xlabel('Frequency, Hz'),ylabel('Magnitude, dB')
title('Magnitude response for Case c design')
axis([0,5000,-80,5]), grid
Computing aliasing characteristic for Case c design
subplot(2,1,2)
k=0;
for l=1:5
for r=1:201
FF(r)=(r-1)*Fy/400;
omega=2*pi*FF(r)/Fy;
W2=exp(-i*n*(omega-2*k*pi)/M);
Halc(r)=sum(hc.*W2);
end
k=k+1;
plot(FF,20*log10(abs(Halc)))
hold on
end
xlabel('Frequency, Hz'),ylabel('Magnitude, dB')
title('Aliasing for Case c design')
axis([0,1000,-80,5]), grid on
hold off
disp('END OF EXAMPLE 3.8')
3.5 Sampling Rate Alteration of Bandpass Signals
The spectrum of bandpass signal is concentrated in the narrow band of frequencies above the zero frequency. The central frequency 
of the signal is generally much larger than the signal bandwidth 
, and the spectrum of the bandpass signal can occupy any band in the range 
.
of the signal is generally much larger than the signal bandwidth , and the spectrum of the bandpass signal can occupy any band in the range .The most popular technique for the sampling-rate alteration of bandpass signals is based on the integer-band decimation and interpolation. This technique can be applied when the signal spectrum is confined to the frequency range
where k is a positive integer. The bandwidth of the signal is restricted to 
, and boundary frequencies of the spectrum are the integer multiples of 
.
, and boundary frequencies of the spectrum are the integer multiples of .3.5.1 Bandpass Signal Decimation
In decimation, prior to downsampling , the signal components outside the desired frequency range have to be eliminated.
The signal is filtered first with an antialiasing filter 
as shown in the block diagram. The ideal magnitude response for the integer-band decimation filter is defined by
as shown in the block diagram. The ideal magnitude response for the integer-band decimation filter is defined by
The MATLAB program illustrates the decimation process of the bandpass signal including the antializaing bandpass filter design.
Example 3.9
clear all, close all
Fx=10000; Fy=2000; %Input and output frequencies
M=5; % Decimation factor
Generating the input signal x
F=[0,2100/5000,2800/5000,2900/5000,3100/5000,3400/5000,3800/5000,...
4000/5000,4500/5000,1];
A=[0,0.05,1,0,0,0.7,0.8,0,0.6,0];
x=fir2(511,F,A); % Generating the input signal signal x[n]
[X,f]=freqz(x,1,1024,10000);
figure
subplot(4,1,1)
plot(f,abs(X)), ylabel('|X(e^j^\omega)|')
Bandpass filter design. Filter is specified to extract the frequency band 2000-3000 Hz (k=2)
hBP=firgr(68,[0,1800/5000,2200/5000,2800/5000,3200/5000,1],[0,0,1,1,0,0]);
[HBP,f]=freqz(hBP,1,1024,10000); % Bandpass filter frequency response
subplot(4,1,2), plot(f,20*log10(abs(HBP))), ylabel('|H_B_P(e^j^\omega)|')
axis([0,5000,-60,5])
Bandpass filtering to obtain intermediate signal v
v=filter(hBP,1,x); % Bandpass filtering to obtain intermediate signal v[n]
[V,f]=freqz(v,1,1024,10000); % Spectrum of the intermediate signal v[n]
subplot(4,1,3), plot(f,abs(V)), ylabel('|V(e^j^\omega)|')
Factor-of-M downsampling
y=downsample(v,5); % Factor-of-M downsampling
Spectrum of decimated signal y
[Y,ff]=freqz(y,1,1024,2000); % Spectrum of decimated signal
subplot(4,1,4); plot(ff,abs(Y)), ylabel('|Y(e^j^\omega)|')
xlabel('Frequency, Hz')
disp('END OF EXAMPLE 3.9')
3.5.2 Integer Band Interpolation
The integer-band interpolation means extracting th kth image from L spectral images of an up-sampled-by-L signal. The structure of an integer bandpass interpolator consisting of a cascade connection of up-sampler and bandpass filter is shown in the Figure.
Example 3.10
clear all,close all
Fx=2000; Fy=10000; %Input and output frequencies
L=5; % Interpolation factor
Generating the input signal x
F=[0,100/1000,800/1000,900/1000,1]; A=[0,0,1,0,0];
x=fir2(102,F,A); % generating the input signal signal x[n]
figure
[X,ff]=freqz(x,1,1024,2000); %Spectrum of the input signal
subplot(4,1,1); plot(ff,abs(X)), ylabel('|Y(e^j^\omega)|')
xlabel('Frequency, Hz')
Bandpass filter design. The bandpass filter is specified to extract the frequency band 3000- 4000 Hz (k = 3).
hBP=5*firgr(64,[0,2800/5000,3200/5000,3800/5000,4200/5000,1],[0,0,1,1,0,0]);
Factor-of-M upsampling to obtain intermediate signal v
v=upsample(x,5);% upsampling
[V,f]=freqz(v,1,1024,10000); % Spectrum of the intermediate signal v[m]
subplot(4,1,2)
plot(f,abs(V)), ylabel('|V(e^j^\omega)|')
[HBP,f]=freqz(hBP,1,1024,10000); % Frequency response of the bandpass filter
subplot(4,1,3)
plot(f,20*log10(abs(HBP))), ylabel('|H_B_P(e^j^\omega)|')
axis([0,5000,-60,5])
Band pass filtering to obtain interpolated signal y
y=filter(hBP,1,v); % filtering
[Y,f]=freqz(y,1,1024,10000);
subplot(4,1,4)
plot(f,abs(Y)), ylabel('|Y(e^j^\omega)|')
xlabel('Frequency, Hz')
disp('END OF EXAMPLE 3.10')
disp(' END OF CHAPTER III')