Untitled

by Ljiljana Milic

Supplemental material for Chapter V:

5. IIR Filters for Sampling Rate Conversion

Table of Contents

5. IIR Filters for Sampling Rate Conversion Example 5.1 5.2 IIR Filter Structures Based on Polyphase Decomposition Example 5.2 5.3 The Role of Extra Filter Example 5.3 5.4 Polyphase IIR Structure with Two All-Pass Subfilters: IIR Halfband Filter Example 5.4 5.5 IIR structures with Two All-pass Subfilters: Applications of EMQF Filters Example 5.5

Example 5.1

Design of example IIR and FIR filters to satisfy Case a tolerance scheme for the factor-of-5 sampling rate conversion.

Specifications :

Passband and stopband cutoff frequencies:

Passband and stopband ripples:

clear all, close all
ap=0.1; as=60;
Fp=0.09; Fs=1/5;

FIR filter design

dev = [(10^(ap/20)-1)/(10^(ap/20)+1)  10^(-as/20)]; 
F = [Fp 1/5];   % Cutoff frequencies
A = [1 0];      % Desired amplitudes
[Nord,Fo,Ao,W] = firpmord(F,A,dev);
h = firpm(Nord+4,Fo,Ao,W); % Filter coefficients
[Hfir,f] = freqz(h,1,1024,2); % FIR filter frequency response
PhaseFIR = unwrap(angle(Hfir));

IIR filter design, 5th order elliptic filter

ap=0.1; as=60;
Fp=0.09; 
[B,A]=ellip(5,ap,as,Fp); % Computing elliptic filter coefficients
[Hiir,f] = freqz(B,A,1024,2); % IIR filter frequency response
PhaseIIR = unwrap(angle(Hiir));

Displaying the magnitude and phase responses

figure

plot(f,20*log10(abs(Hiir)),f,20*log10(abs(Hfir))), grid

legend('IIR filter', 'FIR filter')

xlabel('\omega/\pi'),ylabel('Magnitude, dB')

axis([0,1,-80,5])

figure

plot(f,PhaseIIR,f,PhaseFIR)

xlabel('\omega/\pi'),ylabel('Phase, rad')

legend('IIR filter', 'FIR filter')

disp('END OF EXAMPLE 5.1')

END OF EXAMPLE 5.1

5.2 IIR Filter Structures Based on Polyphase Decomposition

Remainder

Developing a polyphase structure requires the decomposition of the filter transfer function into a set of M(L) polyphase components. The overall transfer function

of a decimation/ interpolation is represented in the form

, or

Efficient polyphase implementation is obtained when the subfilters

are allpass functions

, i.e,

for decimation,

for interpolation.

Solutions based on all-pass polyphase subfilters are restricted to the Case c tolerance scheme.

Efficient commutative implementation structures for decimation and interpolationare are shown bellow.

Commutative polyphase structure of an IIIR factor-of M decimator with all-pass

polyphase branches

Commutative polyphase structure of an IIIR factor-of L interpolator with all-pass

polyphase branches

Example 5.2

IIR polyphase filter for factor-of-5 sampling-rate conversion designed to satisfy Case c toletance scheme.

Specifications: passband cutoff

, stopbands

and

, Minimal stopband attenuation 60 dB.

We use MATLAB functions fminimax and optimset from Optimization Toolbox for optimization, and function elliip from Signal Processing Toolbox.

clear all, close all

Optimization

v0 = rand(1,11)/5; % Generating initial values for fminimax
options = optimset('MinAbsMax',length(1:1000),'TolFun',1e-6,...
    'TolCon',1e-6,'MaxIter',100);
[v,Err,maxfval,exitflag,output] = fminimax(@ErrAllp,v0,...
    [],[],[],[],[],[],[],options);
Local minimum possible. Constraints satisfied.

fminimax stopped because the predicted change in the objective function
is less than the value of the function tolerance and constraints 
are satisfied to within the value of the constraint tolerance.

<stopping criteria details>

Extracting filter poles from vector v

z01 = v(1); z02 = v(2); z03 = v(3);
z11 = v(4); z12 = v(5);
z21 = v(6); z22 = v(7);
z31 = v(8); z32 = v(9);
z41 = v(10); z42 = v(11);

Polyphase branches

a0 = conv([1,-z01],[1,-z02]); a0 = conv(a0,[1,-z03]);
a1 = conv([1,-z11],[1,-z12]);
a2 = conv([1,-z21],[1,-z22]);
a3 = conv([1,-z31],[1,-z32]);
a4 = conv([1,-z41],[1,-z42]);
A0=upsample(a0,5); A1=upsample(a1,5); A2=upsample(a2,5); A3=upsample(a3,5);...
    A4=upsample(a4,5);
B0=fliplr(A0); B1=[0,fliplr(A1)]; B2=[0,0,fliplr(A2)];...
    B3=[0,0,0,fliplr(A3)]; B4=[0,0,0,0,fliplr(A4)]; 

Computing the phase responses of polyphase branches

[PH0,f] = freqz(B0,A0,1000,2); [PH1,f] = freqz(B1,A1,1000,2);...

[PH2,f] = freqz(B2,A2,1000,2); [PH3,f] = freqz(B3,A3,1000,2);...

[PH4,f] = freqz(B4,A4,1000,2);

figure

plot(f,unwrap(angle(PH0)),f,unwrap(angle(PH1)),f,unwrap(angle(PH2)),...

f,unwrap(angle(PH3)),...

f,unwrap(angle(PH4)))

legend('Polyphase branches')

xlabel('Normalized frequency [\omega/\pi]'), ylabel('Phase [rad]')

Overall IIR filter frequency response

H=PH0+PH1+PH2+PH3+PH4;

Equivalent FIR filter

ap = 0.000015; as = 60; Fp = 0.8/5;
dev = [(10^(ap/20)-1)/(10^(ap/20)+1),10^(-as/20),10^(-as/20),10^(-as/20)]; 
F = [Fp,2/5-Fp,2/5+Fp,4/5-Fp,4/5+Fp,1];   % Cutoff frequencies
A = [1,0,0,0];      % Desired amplitudes
[Nord,Fo,Ao,W] = firpmord(F,A,dev); % Estimating the filter order
b = firpm(126,Fo,Ao,W); % Computing the filter coefficients
[Hfir,f] = freqz(b,1,1000,2); % Computing the filter frequency response

Plots of IIR and FIR filter frequency responses

figure

plot(f,20*log10(abs(H)/5),f,20*log10(abs(Hfir))), ylabel('Magnitude [dB]')

axis([0,1,-80,2])

legend('IIR filter','FIR filter')

figure

plot(f,unwrap(angle(H)),f,unwrap(angle(Hfir)))

legend('IIR filter','FIR filter')

xlabel('Normalized frequency [\omega/\pi]'), ylabel('Phase [rad]')

axis([0,1,-50,0])

Plots of passband details: IIR and FIR

figure

%subplot(2,1,1)

plot(f(1:160),20*log10(abs(H(1:160))/5),f(1:160),20*log10(abs(Hfir(1:160))))

%hold on

%figure (3)

%plot(f(1:160),20*log10(abs(Hfir(1:160)))),grid

legend('IIR filter','FIR filter')

xlabel('Normalized frequency [\omega/\pi]'), ylabel('Magnitude, dB')

Generating the input signal {x[n]}

n = 0:510;
x = cos(2*pi*200*n/10000) + cos(2*pi*400*n/10000) + 0.7*cos(2*pi*3800*n/10000) + randn(size(n));
[X,F1] = freqz(x,1,1024,10000);

IIR filter, structure verification: polyphase downsampling by M=5 and filtering

x0 = x(1:5:length(x)); x1 = [0,x(5:5:length(x))]; x2 = [0,x(4:5:length(x))]; 
x3 = [0,x(3:5:length(x))]; x4 = [0,x(2:5:length(x))];
y0 = filter(fliplr(a0)/5,a0,x0); y1 = filter(fliplr(a1)/5,a1,x1);
y2 = filter(fliplr(a2)/5,a2,x2);
y3 = filter(fliplr(a3)/5,a3,x3); y4 = filter(fliplr(a4)/5,a4,x4);
ydec = y0 + y1 + y2 + y3 + y4; % Decimated signal
[Ydec,F2] = freqz(ydec,1,512,2000);

Spectra of input signal and of decimated signal

figure

subplot(2,1,1), plot(F1,abs(X)), legend('Spectrum: input signal')

ylabel('|X(F)|'), axis([0,5000,0,300])

subplot(2,1,2),

plot(F2,abs(Ydec)), legend('Spectrum: decimated signal')

ylabel('|Y_d_e_c(F)|'), xlabel('Frequency [Hz]')

disp('END OF EXAMPLE 5.2')

END OF EXAMPLE 5.2

5.3 The Role of Extra Filter

The single-stage equivalent

is given by:

and

where the extra filter G(z) is designed to improve the overall system frequency response.

Example 5.3

In this example, we examine the frequency response of the sampling rate conversion system shown in the Figure when the sampling rate conversion factor is 5.

Design specifications:

Conversion factor M = 5 and Case a tolerance scheme.
Passband edge frequency =.
Stopband edge frequency = .
Peak passband riple dB.
Minimal stopband attenuation dB.

Solution:

For decimation (interpolation) filter

we take the polyphase realization from Example 5.2. Our choice for

is the 5th order elliptic filter.

clear all, close all

IIR filter poles

z01=-0.03247627480; z02=-0.4519480048; z03=-0.9477051753;
z11=-0.08029157130; z12=-0.5548998293;
z21=-0.1417079348; z22=-0.6883346404;
z31=-0.2320513100; z32=-0.7961481351;
z41=-0.3532045984; z42=-0.8755417392;

Filter

: polyphase branches

a0=conv([1,-z01],[1,-z02]); a0=conv(a0,[1,-z03]);
a1=conv([1,-z11],[1,-z12]);
a2=conv([1,-z21],[1,-z22]);
a3=conv([1,-z31],[1,-z32]);
a4=conv([1,-z41],[1,-z42]);
A0=upsample(a0,5); A1=upsample(a1,5); A2=upsample(a2,5);
A3=upsample(a3,5); A4=upsample(a4,5);
B0=fliplr(A0); B1=[0,fliplr(A1)]; B2=[0,0,fliplr(A2)];
B3=[0,0,0,fliplr(A3)]; B4=[0,0,0,0,fliplr(A4)]; 

Computing the frequency responses of polyphase branches

[PH0,f]=freqz(B0,A0,1024,2); [PH1,f]=freqz(B1,A1,1024,2);
[PH2,f]=freqz(B2,A2,1024,2); [PH3,f]=freqz(B3,A3,1024,2); 
[PH4,f]=freqz(B4,A4,1024,2);

frequency response

H1=PH0+PH1+PH2+PH3+PH4;

filter design

[b,a]=ellip(5,0.1,60,0.8);
[G,f]=freqz(b,a,1024,2);
bb=upsample(b,5);
aa=upsample(a,5);
[Gu,f]=freqz(bb,aa,1024,2);

Computing frequency response of

H=(H1/5).*Gu;

figure

plot(f,20*log10(abs(G)),'r')

ylabel('Gain, dB'),xlabel('\omega/\pi')

title('Extra ilter')

axis([0,1,-80,5])

grid

figure

subplot(2,1,1)

plot(f,20*log10(abs(H1)/5),'b'), grid

title('Decimation filter H_1(z) and periodic extra filter G(z^5)')

hold on

plot(f,20*log10(abs(Gu)))

ylabel('Gain, dB')

axis([0,1,-80,5])

subplot(2,1,2)

plot(f,20*log10(abs(H)),'k')

title('Overall filter H(z)')

xlabel('\omega/\pi')

axis([0,1,-80,5])

grid

figure

subplot(3,1,1)

plot(f(1:104),20*log10(abs(H(1:104))),'k'), grid

xlabel('\omega/\pi'),ylabel('Gain, dB')

title('Passband characteristic')

axis([0,0.1,-0.1,0.01])

disp('END OF EXAMPLE 5.3')

END OF EXAMPLE 5.3

5.4 Polyphase IIR Structure with Two All-Pass Subfilters:

IIR Halfband Filter

The transfer function of a lowpass IIR halfband filter

can be represented as a sum of two all-pass subfilters

Here,

are polyphase components of the transfer function

Efficient polyphase structures for decimation and interpolation for M = L = 2 are shown bellow:

(a) Polyphase factor-of-2 decimator. (b) Polyphase factor-of-2 interpolator.

Example 5.4

close all, clear all

Design and polyphase representation of halfband filter

p = [0.4863i,-0.4863i,0.8453i,-0.8453i];  % Halfband filter poles
% Coefficients of allpass branches
beta0 = (abs(p(1)))^2;    beta1 = (abs(p(3)))^2; 

All-pass branches

[A0,f] = freqz([beta0,0,1],[1,0,beta0],512,2); 
[A1,f] = freqz([0,beta1,0,1],[1,0,beta1],512,2);

Half-band filter frequency response

H = (A0 + A1)/2; % Computing the IIR halfband filter frequency response

figure

plot(f,20*log10(abs(H))), grid

xlabel('\omega/\pi'),ylabel('Gain, dB')

axis([0,1,-60,2])

text(0.65,-20,'Filter Characteristic')

Generating the input signal

F = [0,0.1,0.46,1]; A=[0,1,0,0]; x1 = fir2(256,F,A);

x = x1 + 0.008*cos(2*pi*0.35*(0:256));% Original signal

X = fft(x); % Spectrum of the original signal

figure

plot((0:127)/128,abs(X(1:128))),grid

xlabel('\omega/\pi'),ylabel('|X|')

text(0.40,0.8,'Original signal')

axis([0,1,0,1.2])

Decimation-by-2

u0 = x(1:2:length(x));

u1 = [0,x(2:2:length(x)-1)];

y0 = filter([beta0,1],[1,beta0],u0);

y1 = filter([beta1,1],[1,beta1],u1);

ydec = (y0 + y1)/2;

Ydec = fft(ydec);

figure

plot((0:63)/64,abs(Ydec(1:64))),grid

xlabel('\omega/\pi'),ylabel('|Y_{dec}|')

text(0.65,0.4,'Decimated signal')

axis([0,1,0,0.6])

Interpolation-by-2

y0 = filter([beta0,1],[1,beta0],ydec);

y1 = filter([beta1,1],[1,beta1],ydec);

yy = [y0;y1];

yint = yy(:);

Yint = fft(yint);

figure

plot((0:127)/128,abs(Yint(1:128))),grid

text(0.65,0.8,'Interpolated signal')

xlabel('\omega/\pi'), ylabel('|Y_{int}|'), axis([0,1,0,1.2])

disp('END OF EXAMPLE 5.4')

END OF EXAMPLE 5.4

5.5 IIR structures with Two All-pass Subfilters:

Applications of EMQF Filters

Example 5.5

Remainder

Consider two solutions for factor-of-8 decimator implemented in three stages as indicated in the Figure.

(a) Filters

, and

are 5th order halfband filters.

(b) Filters

are 5th order halfband filters, and

is an EMQF (Elliptic Minimum Q-Factor Filter) designed to achieve Case a specifications for M = 2.

Three-stage decimator for M = 8

The single-stage equivalent of the overall transfer function of the decimator is given by

close all, clear all

Halfband filter design

p=[0.4863i,-0.4863i,0.8453i,-0.8453i];  % Halfband filter poles
beta2HB=(abs(p(1)))^2;    beta3HB=(abs(p(3)))^2; % Coefficients
[A0HB,w]=freqz([beta2HB,0,1],[1,0,beta2HB],1000); % Branch A0
[A1HB,w]=freqz([0,beta3HB,0,1],[1,0,beta3HB],1000); % Branch A1
[A0HBu,w]=freqz(upsample([beta2HB,0,1],2),upsample([1,0,beta2HB],2),1000); 
[A1HBu,w]=freqz(upsample([0,beta3HB,0,1],2),upsample([1,0,beta3HB],2),1000);
H_HBu=(A0HBu+A1HBu)/2; 
[A0HBu4,w]=freqz(upsample([beta2HB,0,1],4),upsample([1,0,beta2HB],4),1000); 
[A1HBu4,w]=freqz(upsample([0,beta3HB,0,1],4),upsample([1,0,beta3HB],4),1000);
H_HBu4=(A0HBu4+A1HBu4)/2; 
H_HB=(A0HB+A1HB)/2; % Half-band filter frequency response
figure 
plot(w/pi,20*log10(abs(H_HB)))
hold on

EMQF filter design

N=7;
alpha=-(1/8+1/32);
f3db=acos(-alpha)/(2*pi);
ksi=(tan(2*pi*0.25/2)/tan(2*pi*f3db/2))^2;
omegap=2*atan(1/sqrt(ksi));
fp=omegap/(2*pi);
[b,a,z,p,k] = halfbandiir(7,2*fp);
b=sort(abs(p));
betaHB=b(2:2:N).^2;
alpha=-cos(2*pi*f3db);
alpha1=(1-sqrt(1-alpha^2))/alpha;
beta=(betaHB+alpha1^2)./(betaHB*alpha1^2+1);
a1=alpha1;
p1 = [1 a1];
for ind = 2:2:(N-1)/2
  p1 = conv(p1,[1 alpha*(1+beta(ind)) beta(ind)]);
end
p2 = 1;
for ind = 1:2:(N-1)/2
  p2 = conv(p2,[1 alpha*(1+beta(ind)) beta(ind)]);
end

Allpass branches

Hab=p1;
Hai=fliplr(Hab);
Hbb=p2;
Hbi=fliplr(Hbb);
[Ha,f]=freqz(Hab,Hai,1000,2);
[Hb,f]=freqz(Hbb,Hbi,1000,2);

EMQF filter frequency response

H_EMQF=(Ha+Hb)/2;

plot(f,20*log10(abs(H_EMQF)),'r')

legend('5th ord. HB filter','7th ord. EMQF filter')

legend('boxoff')

ylabel('Gain, dB'), axis([0,1,-60,5]), grid

Hab_u=upsample(Hab,4);

Hai_u=upsample(Hai,4);

Hbb_u=upsample(Hbb,4);

Hbi_u=upsample(Hbi,4);

[Ha_u,f]=freqz(Hab_u,Hai_u,1000,2);

[Hb_u,f]=freqz(Hbb_u,Hbi_u,1000,2);

H_EMQF_u=(Ha_u+Hb_u)/2;

Single-stage equivalent for solution (a)

H_SINGLE_HB=H_HB.*H_HBu.*H_HBu4; % Solution (a)

Single-stage equivalent for solution (b)

H_SINGLE=H_HB.*H_HBu.*H_EMQF_u; % Solution (b)

Displaying gain responses for solutiona (a) and (b)

figure

plot(f,20*log10(abs(H_SINGLE_HB)),f,20*log10(abs(H_SINGLE)))

grid

legend('Solution (a)','Solution (b)')

legend('boxoff')

xlabel('\omega/\pi'), ylabel('Gain, dB')

axis([0,1,-80,5])

g=axes('Position',[0.35 0.60 0.40 0.12]);

plot(f(1:102),20*log10(abs(H_SINGLE_HB(1:102))),...

f(1:102),20*log10(abs(H_SINGLE(1:102))))

grid

axis([0,0.1,-3*0.001,0.5*0.001])

disp('END OF EXAMPLE 5.5')

END OF EXAMPLE 5.5

disp( ' END OF CHAPTER V')

                      END OF CHAPTER V

function Err = ErrAllp(v)
% Extracting poles from vector v
z01 = v(1); z02 = v(2); z03 = v(3);
z11 = v(4); z12 = v(5);
z21 = v(6); z22 = v(7);
z31 = v(8); z32 = v(9);
z41 = v(10); z42 = v(11);
% Polyphase branches
a0 = conv([1,-z01],[1,-z02]); a0 = conv(a0,[1,-z03]);
a1 = conv([1,-z11],[1,-z12]);
a2 = conv([1,-z21],[1,-z22]);
a3 = conv([1,-z31],[1,-z32]);
a4 = conv([1,-z41],[1,-z42]);
A0=upsample(a0,5); A1=upsample(a1,5); A2=upsample(a2,5);...
    A3=upsample(a3,5); A4=upsample(a4,5);
B0=fliplr(A0); B1=[0,fliplr(A1)]; B2=[0,0,fliplr(A2)];...
    B3=[0,0,0,fliplr(A3)]; B4=[0,0,0,0,fliplr(A4)];
[PH0,f] = freqz(B0,A0,1000,2); [PH1,f] = freqz(B1,A1,1000,2);...
    [PH2,f] = freqz(B2,A2,1000,2); [PH3,f] = freqz(B3,A3,1000,2);...
    [PH4,f] = freqz(B4,A4,1000,2);
H = PH0+PH1+PH2+PH3+PH4;
H = abs(H)/5;
Err = [(1-H(1:160)'),zeros(size(161:240)),(H(241:560))',...
    zeros(size(561:640)),(H(641:960))',zeros(size(961:1000))];
end
%--------------------------------------------------------------------EOF