%
% 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 \left\lbrack 1-\delta_p ,1\right\rbrack.
%
%% 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 H\left(z^L \right), and the second filter G\left(z\right) 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
% H_{\textrm{IFIR}} \left(z\right)=H\left(z^L \right)G\left(z\right)
%
% Design IFIR filter to satisfy Case a tolerance scheme.
% Design specifications:
% Sampling-rate conversion factor M = 5.
% Passband edge \omega_p =0\ldotp 09\pi. Stopband edge at \omega_s =\frac{\pi }{5}.
% Pasband and stopband ripples are given in decibels with a_p =0\ldotp 1\;\textrm{dB}, and a_s =60\;\textrm{dB}.

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 H\left(z^L \right) and \textrm{image}\;\textrm{suppression}\;\textrm{filter}\;G\left(z\right)

[h,g] = ifir(L,'low',[Fp, Fs],dev); 

% Computing frequency responses of H\left(z^L \right) and G\left(z\right)

[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 H\left(z\right) 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 Y\left(e^{j\omega } \right) is expressible in terms of the spectrum of the high-rate signal X\left(\left(e^{j\omega } \right.\right) and the frequency response of the antialiasing filter H\left(e^{j\omega } \right)
% Y\left(e^{j\omega } \right)=\frac{1}{M}\sum_{k=0}^{M-1} X\left(e^{j\left(\omega -2\pi k\right)/M} \right)H_{\textrm{al}}^k \left(e^{j\omega } \right)
% where
% H_{\textrm{al}}^k \left(e^{j\omega } \right)=H\left(e^{j\left(\omega -2\pi k\right)/M} \right),\;\;\;k=0,1,\cdots ,M-1
% Here, for k = 0 function H_{\textrm{al}}^k \left(e^{j\omega } \right) represents the unalized characteristic, whereas H_{\textrm{al}}^k \left(e^{j\omega } \right) 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 H_{\textrm{al}}^k \left(e^{j\omega } \right) and usually express in decibels,
% a_{\textrm{al}}^k \left(\omega \right)=20\log_{10} \left|H_{\textrm{al}}^k \left(e^{j\omega } \right)\right|,\;\;\;k=0,1,\cdots ,M-1
% 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 \omega_c of the signal is generally much larger than the signal bandwidth \Delta_B, and the spectrum of the bandpass signal can occupy any band in the range \left\lbrack 0,\pi \right\rbrack.
% 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
% \frac{k\pi }{M}\le \omega \le \frac{\left(k+1\right)\pi }{M}
% where k is a positive integer. The bandwidth of the signal is restricted to \Delta_B \le \pi /M, and boundary frequencies of the spectrum are the integer multiples of \pi /M.
%
%% 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 H_{\textrm{BP}} \left(z\right)\;as shown in the block diagram.  The  ideal magnitude response for the integer-band decimation filter is defined by
% \left|H_{\textrm{BP}} \left(e^{j\omega } \right)\right|=\left\lbrace \begin{array}{ll}
% 1, & k\pi /M\le \left|\omega \right|\le \left(k+1\right)\pi /M\\
% 0, & \textrm{otherwise}
% \end{array}\right.
% 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')