%%
% by Ljiljana Milic
% Supplemental matherial 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 H_K \left(z\right)

hk = firgr(50,[0,2*Fp/(F0/2),2*Fs/(F0/2),1],[1,1,0,0]);  

% Designing the decimation filter H_D \left(z\right)

hd = firhalfband('minorder',(2*Fs + Fp)/(3*F0/2),delta); 

% Interpolation filter H_I \left(z\right)

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 H_D \left(z\right)\;and H_I \left(z\right), and an eliptic mimimal Q-factors (EMQF) filter for the kernel filter H_K \left(z\right).

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 H_K \left(z\right)
% 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 H_D \left(z\right)
% 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:
% \left(F_s^2 -F_p^2 \right)M^3 -{\left(F_s +F_p \right)}^2 M^2 +2F_0 {\left(F_s +F_p \right)}^M -F_0^2 =0
%% Example 9.3
% Specifications for FIR lowpass filter:
% Bandpass edge frequency: F_p =50\;\textrm{Hz}
% Stopband edge frequency: F_s =100\;\textrm{Hz}
% Sampling frequency: F_0 =2000\;\textrm{Hz}
% Ripple tolerance in the passband: \delta_1 =0\ldotp 01\;
% Minimal stopband attenuation 60 dB (\delta_2 =0\ldotp 001)

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 H_K \left(z\right), decimation filter H_I \left(z\right), interpolation filter H_I \left(z\right), and delay element z^{-D}.
%
% Wideband filter specifications:
% Sanpling frwquency: F_0 =2000\;\textrm{Hz}
% Passband edge frequency F_p =700\;\textrm{Hz}
% Stopband edge frequency F_s =770\;\textrm{Hz}
% Minimal stopband attenuation a_s =80\;\textrm{dB}
% 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: F_0 =4000\;\textrm{Hz}
% Passband edge frequency F_p =700\;\textrm{Hz}
% Stopband edge frequency F_s =770\;\textrm{Hz}
% Minimal stopband attenuation a_s =60\;\textrm{dB}
%
%
% 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 H_{\textrm{D1}} \left(z\right)\; and H_{\textrm{I1}} \left(z\right)\;.

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