%%
% by Ljiljana Milic
% Supplemental material for Chapter IV
%% 4. FIR Filters for Sampling Rate Conversion
% Table of Contents
% 4. FIR Filters for Sampling Rate Conversion
%     4.1 Introduction
%     4.2 Structure Verification and Simulation Using MATLAB
%         Example 4.1
%         Example 4.2
%     4.3 Memory Saving Structures for FIR Polyphase Decimators and Interpolators
%         Example 4.3
%     4.4 Computational Efficiency of Single-Stage and Multistage Systems
%         Example 4.4
%% 4.1 Introduction
% The transfer function of an FIR filter is a the polynomial in z^{-1}, and is usually writen in the form
% H\left(z\right)=\sum_{k=0}^{N-1} h\left\lbrack k\right\rbrack z^{-k}
% Here, we express H\left(z\right) by means of M (or L) polyphase branches E_0 \left(z\right),E_1 \left(z\right),\cdots {,E}_M \left(z\right)
% \;H\left(z\right)=\sum_{k=0}^{M-1} z^{-1} E_k \left(z^M \right)
% where
% E_k \left(z\right)=\sum_n^{\left\lfloor N/M\right\rfloor } h\left\lbrack \textrm{Mn}+k\right\rbrack z^{-n}
% Based on the polyphase implementation, operations in sampling-rate conversion systems can be executed at the lower sampling rate. The efficient commutative implementation structures for polyphase decimator| and interpolator are shown and verified in the next subsections.
%
%% 4.2 Structure Verification and Simulation Using MATLAB
%
%
% Commutative polyphase structure of a decimator.
%
%
% Commutative polyphase structure of an interpolator.
%% Example 4.1
% Decimation and interpolation by integer factors M = L = 5
% Design of optimal linear-phase FIR filter
% Creation of the input signal
% Polyphase decomposition to 5 polyphase components
% Polyphase decimation according to the Figure shown above
% Polyphase interpolation according to the Figure shown above
% Structure verification

clear all, close all

% STEP 1: Lowpass FIR filter design to satisfy Case a tolerance scheme.
% Filter specifications: Fp = 0.08, Fs = 1/M - Fp, Passband ripple ap = 0.01 dB, minimum stopband
%  attenuation as = 60 dB.Specifications are met with an optimal FIR filter of the length N = 60.

h = firpm(59,[0,2*400/10000,1000/5000,1],[1,1,0,0]); % Filter coefficients
[H,F1] = freqz(h,1,1024,10000);  % Computation of filter frequency response
figure 
subplot(4,1,1), plot(F1,20*log10(abs(H))), legend('Filter magnitude response'),...
    ylabel('Gain, dB'),axis([0,5000,-80,5])

% STEP 2: Generating the input signal
% The input signal is a causal sequence formed as a sum of three sinusoidal sequencies corrupted in noise.
% x\left\lbrack n\right\rbrack =\cos \left(2\pi F_1 n/F_x \right)+\cos \left(2\pi F_2 n/F_x +0\ldotp 7\cos \left(2\pi F_3 n/F_x \right)+s\left\lbrack n\right\rbrack \right)
% where F_1 =200\;\textrm{Hz}, F_2 =400\;\textrm{Hz}, F_3 =3800\;\textrm{Hz}, s\left\lbrack n\right\rbrack is random signal with the normal distribution.
% The sampling frequency at the input F_x =10000\;\textrm{Hz}. The sequence length N = 511 samples.

n = 0:510; % Time index
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); % Computing the signal spectrum
subplot(4,1,2), plot(F1,abs(X)),... 
legend('Spectrum: input signal'), ylabel('|X(F)|'), axis([0,5000,0,300])

% STEP 3: Polyphase decomposition

E = reshape(h,5,length(h)/5);

% STEP 4: Polyphase downsampling 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(E(1,:),1,x0); y1 = filter(E(2,:),1,x1); y2 = filter(E(3,:),1,x2);...
    y3=filter(E(4,:),1,x3); y4=filter(E(5,:),1,x4);
ydec = y0 + y1 + y2 + y3 + y4; % Decimated signal
[Ydec,F2] = freqz(ydec,1,512,2000);
subplot(4,1,3), plot(F2,abs(Ydec)),legend('Spectrum: decimated signal'),...
     ylabel('|Y_d_e_c(F)|'), axis([0,1000,0,60])

% STEP 5: Poliphase filtering and up-sampling

x = ydec; % Input to the interpolator
u0 = 5*filter(E(1,:),1,x); u1 = 5*filter(E(2,:),1,x); u2 = 5*filter(E(3,:),1,x);
u3 = 5*filter(E(4,:),1,x); u4 = 5*filter(E(5,:),1,x);
U = [u0;u1;u2;u3;u4];
yint = U(:); % interpolated signal
[Yint,F3] = freqz(yint,1,1024,10000);
subplot(4,1,4), plot(F3,abs(Yint)), legend('Spectrum: interpolated signal'),...
    xlabel('Frequency, Hz'), ylabel('|Y_i_n_t(F)|')
axis([0,5000,0,300])
disp('END OF EXAMPLE 4.1')

%
%% Example 4.2
% Decimation and interpolation by integer factors M = L = 5
% Design of optimal linear-phase FIR filter
% Creation of the input signal
% Polyphase decomposition to 5 polyphase components
% Simulation of real-time decimation
% Simulation of real-time interpolation.

clear all, close all

% STEP 1: Lowpass FIR filter design  to satisfy Case a tolerance scheme.
% Filter specifications from example 4.1.

h = firpm(59,[0,400/5000,1000/5000,1],[1,1,0,0]); % Filter coefficients
[H,F1] = freqz(h,1,1024,10000); % Computing the signal spectrum
figure 
subplot(4,1,1), plot(F1,20*log10(abs(H))),...
    legend('Filter magnitude response'), ylabel('Gain, dB'),axis([0,5000,-80,5])

% STEP 2: Generating the input signal as given in Example 4.1

n = 0:510;  % Time index
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); % Computing the spectrum of the input signal
subplot(4,1,2)
plot(F1,abs(X)), legend('Spectrum: input signal'), ylabel('|X(F)|')
axis([0,5000,0,300])

% STEP 3
% Polyphase decomposition

% Polyphase decomposition

E = reshape(h,5,length(h)/5);

% STEP 4: Simulation of decimator

vk_i = zeros(size(E)); %Initial states of polyphase components
rot = 0; % Initial rotator position
ydec = [];
for n=1:length(x)
xn = x(n);
vk = vk_i(rot+1,:); ek = E(rot+1,:);
vk = [xn,vk(1:length(vk)-1)]; yn = sum(vk.*ek); % filtering
vk_i(rot+1,:)=vk;
if n==1, yy=0; else, end
if rot == 0 ss = 4; yy = yy+yn; ydec = [ydec,yy]; else end
if rot == 1 ss = 0; yy = yy+yn; else end
if rot == 2 ss = 1; yy = yy+yn; else end
if rot == 3 ss = 2; yy = yy+yn; else end
if rot == 4 ss = 3; yy = yn; else end
rot = ss;
end
[Ydec,F2] = freqz(ydec,1,512,2000);
subplot(4,1,3), plot(F2,abs(Ydec)),...
    legend('Spectrum: decimated signal'), ylabel('|Y_d_e_c(F)|')
axis([0,1000,0,60])

% STEP 5: Simulation of interpolator

yi = [0]; xk_i = zeros(size(E));
rot = 0;
k = 1;
for n=1:length(ydec)*5
    xn=ydec(k);
xk = xk_i(rot+1,:); ek = E(rot+1,:); 
xk=[xn,xk(1:length(xk)-1)]; yn=sum(xk.*ek); % filtering
xk_i(rot+1,:)=xk; yn=5*yn;
if rot == 0 ss = 1; yi = [yi,yn]; else end
if rot == 1 ss = 2; yi = [yi,yn]; else end
if rot == 2 ss = 3; yi = [yi,yn]; else end
if rot == 3 ss = 4; yi = [yi,yn]; else end
if rot == 4 ss = 0; k = k+1; yi = [yi,yn]; else end
rot = ss;
end
[Yint,F3] = freqz(yi,1,1024,10000);
subplot(4,1,4)
plot(F3,abs(Yint)); legend('Spectrum: interpolated signal')
xlabel('Frequency, Hz'),ylabel('|Y_i_n_t(F)|')
axis([0,5000,0,300])
disp('END OF EXAMPLE 4.2')

%
%% 4.3 Memory Saving Structures for FIR Polyphase Decimators and Interpolators
%
% Memory saving structures based on polyphase implementation of decimator and interpolator are shown bellow for sampling-rate conversion factor M = L = 3.
%
% Memory saving structure of the polyphase factor-of-3 decimator.
%
% Memory saving structurefor the polyphase factor-of-3 interpolator.
%
%% Example 4.3
% This exampl  is developed to simulate the memory saving commutative decimator and interpolator for the sampling-rate conversion factor M = L = 3. The sampling-frequency of the test signal at the input amounts 10000 Hz.
% Program is organized in 5 steps.

clear all, close all

% STEP 1:  Lowpass FIR filter design:
% Optimal FIR filter of the length N = 12. Specifications: Case b tolerance scheme.

h = firpm(11,[0,800/5000,2*1666/5000-800/5000,1],[1,1,0,0]);
[H,F1] = freqz(h,1,1024,10000);
figure 
subplot(4,1,1), plot(F1,20*log10(abs(H))), legend('Filter magnitude response'), ylabel('Gain, dB')
axis([0,5000,-60,5])

% STEP 2:  Polyphase decomposition

E = reshape(h,3,length(h)/3); E = fliplr(E);

% STEP 3: Input signal

n = 0:510; % Time index
x = cos(2*pi*350*n/10000) + cos(2*pi*750*n/10000) + ...
    0.7*cos(2*pi*3800*n/10000) + randn(size(n));
[X,F1] = freqz(x,1,1024,10000); % Computing the spectrum of the input signal
subplot(4,1,2)
plot(F1,abs(X)), legend('Spectrum: input signal'), ylabel('|X(F)|')
axis([0,5000,0,300])

% STEP 4
% Setting the initial states for decimator

xk = zeros(size(1:length(h)/3));

% Decimation

rot = 0; % Initial rotator position
ydec = [];
for n=1:length(x)
xn = x(n); xk = xn*E((rot+1),:)+xk;
if rot == 0 ss = 2; yn = xk(length(xk)); ydec = [yn,ydec];...
        xk = [0,xk(1:length(xk)-1)]; else end
if rot == 1 ss = 0; else end
if rot == 2 ss = 1; else end
rot = ss; end
[Ydec,F2] = freqz(ydec,1,512,3332);
subplot(4,1,3), plot(F2,abs(Ydec)), legend('Spectrum: decimated signal'),...
    ylabel('|Y_{dec}(F)|')
axis([0,1662,0,100])

% STEP 5: Interpolation

E = 3*fliplr(E);
yi = [];
xk = zeros(size(1:length(h)/3));;
rot = 0; % Initial rotator position
k = 1;
for n=1:length(ydec)*3
    xn = ydec(k); ek = E(rot+1,:);
    if k == 1,  xk = [xn,xk(1:length(xk)-1)]; else end
    yn = sum(xk.*ek); yi = [yi,yn];
if rot == 0 ss = 1; else end                                                                                                                                                                                                                                                                        
if rot == 1 ss = 2; else end
if rot == 2 ss = 0; k = k+1; xk = [xn,xk(1:length(xk)-1)]; else end
rot = ss;
end
[Yint,F3] = freqz(yi,1,1024,10000);
subplot(4,1,4)
plot(F3,abs(Yint)); legend('Spectrum: interpolated signal')
xlabel('Frequency, Hz'),ylabel('|Y_{int}(F)|')
axis([0,5000,0,300])
disp('END OF EXAMPLE 4.3')

% \times
%
%% 4.4 Computational Efficiency of Single-Stage and Multistage Systems
%% Example 4.4
% Implementation of 2-stage factor-of-15 decimator where M=M_1 \times M_2, M_1 =5, M_2 =3.
% Sampling rate conversion from 30 kHz to 2 kHz.
%
% Two-stage decimator: (a) Two atage implementation for M=5\times 3,
% (b)  Single-stage equivalent for M=5\times 3 (developed ftom the 3rd identity).

close all, clear all

% Filter H_1 \left(z\right)

h1 = firpm(18,[0,0.5/15,5/15,7/15,11/15,13/15,1,1],[1,1,0,0,0,0,0,0],[1,5,5,5]);
[H1,F1] = freqz(h1,1,1024,30); % Frequency response  of H1(z)

% Filter H_2 \left(z\right)

h2 = firpm(35,[0,5*0.5/15,5*1/15,1],[1,1,0,0],[1,5]); % Filter H2(z)
[H2,F2] = freqz(h2,1,1024,6); % Frequency response of H2(z)

% Periodic filter H_2 \left(z^5 \right)

h2u = upsample(h2,5); % Periodic filter H_2(z^5)
[H2u,F1] = freqz(h2u,1,1024,30); % Spectrum of the periodic filter H_2(z^5)

% The overall frequency response H\left(z\right)

H = H1.*H2u; % The overall frequency response
figure 
subplot(4,1,1), plot(F2,20*log10(abs(H2))), ylabel('|H_2(F)|'), axis([0,3,-80,2])
title('Magnitude responses, dB')
subplot(4,1,2), plot(F1,20*log10(abs(H2u))), ylabel('|H_2_u(F)|'),axis([0,15,-80,2])
subplot(4,1,3), plot(F1,20*log10(abs(H1))), ylabel('|H_1(F)|'), axis([0,15,-80,2])
subplot(4,1,4), plot(F1,20*log10(abs(H))), ylabel('|H(F)|'), axis([0,15,-80,2])
xlabel('Frequency, kHz')

% Creating the input signal:
%

ff = [0,1/15,1.5/15,0.6,0.8,1]; a = [0,0,0.8,0.7,0.9,0];
x = fir2(1000,ff,a);
x = 5*x+0.01*sin(2*pi*(1:1001)*300/30000) + 0.01*sin(2*pi*(1:1001)*400/30000);
[X,F1] = freqz(x,1,2048,30); % Spectrum of the input signal
figure 
subplot(3,1,1), plot(F1,abs(X)), ylabel('|X(F)|'), legend('Original')
title('Signal spectrum: original, 1st stage, 2nd stage')

%  Polyphase decomposition M_1 =5

h1 = [h1,0];
E1 = reshape(h1,5,length(h1)/5);

% Polyphase  downsampling and filtering, M_1 =5

v0 = x(1:5:length(x)); v1 = [0,x(5:5:length(x))]; v2 = [0,x(4:5:length(x))]; 
v3 = [0,x(3:5:length(x))]; v4 = [0,x(2:5:length(x))];
y0 = filter(E1(1,:),1,v0); y1 = filter(E1(2,:),1,v1); y2 = filter(E1(3,:),1,v2);
y3 = filter(E1(4,:),1,v3); y4 = filter(E1(5,:),1,v4);
yd1 = y0 + y1 + y2 + y3 + y4; % Decimated signal
[Yd1,F2]=freqz(yd1,1,2048,6);
subplot(3,1,2), plot(F2,abs(Yd1)),
legend('Decimated, 1st stage'), ylabel('|Y_d_1(F)|')

%  Polyphase decomposition M_2 =3

E2 = reshape(h2,3,length(h2)/3);

% Polyphase  downsampling and filtering, M_2 =3

x = [yd1,0];
v0 = x(1:3:length(x)); v1 = [0,x(3:3:length(x))]; v2 = [0,x(2:3:length(x))];
y0 = filter(E2(1,:),1,v0); y1 = filter(E2(2,:),1,v1);
y2 = filter(E2(3,:),1,v2);
y = y0 + y1 + y2; % decimated signal
[Y,F3] = freqz(y,1,2048,2);
subplot(3,1,3), plot(F3,abs(Y)), legend('Decimated, 2nd stage'), ylabel('|Y(F)|')
xlabel('Frequency, kHz')

% Comment: Efficiency of implementation
% Single-stage factor-of-15 FIR decimator needs 180 multiplication constanta
% Two-stage 5x3 FIR decimator needs 55 multiplication constants as shown above.

disp('END OF EXAMPLE 4.4')
disp('                    END OF CHAPTER IV')

%