%%
% 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: F_p =0\ldotp 09, F_s =1/5
% Passband and stopband ripples:a_p =0\ldotp 1\;\textrm{dB}, a_s =60\;\textrm{dB}

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')

%
%% 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 H\left(z\right) of a decimation/ interpolation  is represented in the form H\left(z\right)=\frac{1}{M}\sum_{k=0}^{M-1} z^{-k} H_k \left(z^M \right), or   H\left(z\right)=L\sum_{k=0}^{L-1} z^{-k} H_k \left(z^L \right)
% Efficient polyphase implementation is obtained when the subfilters H_k \left(z\right)are allpass functionsA_k \left(z\right), i.e,
% H\left(z\right)=\frac{1}{M}\sum_{k=0}^{M-1} z^{-k} A_k \left(z^M \right) for decimation,
% H\left(z\right)=\sum_{k=0}^{M-1} z^{-k} A_k \left(z^M \right) 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 \omega_p =0\ldotp 8\pi /5, stopbands \left\lbrace 1\ldotp 2\pi /5-2\ldotp 8\pi /5\right\rbraceand \left\lbrace 3\ldotp 2\pi /5-4\ldotp 8\pi /5\right\rbrace, 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);

% 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')

%
%% 5.3 The Role of Extra Filter
%
%
% The single-stage equivalent H\left(z\right)is given by:
% H\left(z\right)=H_1 \left(z\right)G\left(z^M \right) and H\left(z\right)=H_1 \left(z\right)G\left(z^L \right)
% 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 \omega_{p\;}=\;0\ldotp 8\pi /5.
% Stopband edge frequency \omega_{s\;}= \pi /5.
% Peak passband riple a_p =0\ldotp 1\;dB.
% Minimal stopband attenuation a_s =60\;dB.
% Solution:
% For decimation (interpolation) filter H_1 \left(z\right) we take the polyphase realization from Example 5.2. Our choice for G\left(z\right) 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 H_1 \left(z\right): 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);

% H_1 \left(z\right) frequency response

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

% G\left(z\right) 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\left(z\right)

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')

%
%% 5.4 Polyphase IIR Structure with Two All-Pass Subfilters:
%% IIR Halfband Filter
%
% The transfer function of a lowpass IIR halfband filter H^{\textrm{HB}} \left(z\right) can be represented as a sum of two all-pass subfilters A_0^{\textrm{HB}} \left(z^2 \right)\;\textrm{and}{\;A}_1^{\textrm{HB}} \left(z^2 \right)
% H^{\textrm{HB}} \left(z\right)=\frac{1}{2}\left(A_0^{\textrm{HB}} \left(z^2 \right)+A_1^{\textrm{HB}} \left(z^2 \right)\right)
% Here, A_0^{\textrm{HB}} \left(z^2 \right)\;\textrm{and}{\;A}_1^{\textrm{HB}} \left(z^2 \right) are polyphase components of the transfer function H^{\textrm{HB}} \left(z\right).
% 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 H^{\textrm{HB}} \left(z\right)

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')

%
%% 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 H_1 \left(z\right), H_2 \left(z\right), and H_3 \left(z\right) are 5th order halfband filters.
% (b) Filters H_1 \left(z\right), H_2 \left(z\right) are 5th order halfband filters, and H_3 \left(z\right) 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
% H\left(z\right)=H_1 \left(z\right)H_2 \left(z^2 \right)H_3 \left(z^4 \right)
%
%

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')
disp( '                      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
