%
% by Ljiljana Milic
% Supplemental material for Chapter II:
%% 2.  Basics of Multirate Systems
% Table of Contents
% 2.  Basics of Multirate Systems
%     2.1 Frequency-Domain Characterisation of Down-Sampling
%         Example 2.1
%     2.2 Frequency-Domain Characterisation of Up-Sampling
%         Example 2.2
%     2.3 Spectrum of Decimated Signal
%         Example 2.3
%     2.4 Spectrum of Interpolated Signal
%         Example 2.4
%     2.5 The Six Identities
%         First identity
%         Second identity
%         Third Indentity
%         Example 2.5
%         Fourth Identity
%         Fifth Identity
%         Sixth Identity
%         Example 2.6
%     2.6 Cascading Sampling-Rate Alteration Devices
%         Example 2.7
%     2.7 Sampling-Rate Conversion with the Phase Offset
%         Example 2.8
%     2. 8 Polyphase Decomposition and Reconstruction
%         Example 2.9
%     2.9 Multistage Systems
%         Example 2.10
%% 2.1 Frequency-Domain Characterisation of Down-Sampling
%
% Remainder
%
%
% Factor-of-M downsampler
%
% Downsampler's response \left\lbrace y\left\lbrack m\right\rbrack \right\rbraceto the input sequence \left\lbrace x\left\lbrack n\right\rbrack \right\rbrace:
% y\left\lbrack m\right\rbrack =x\left\lbrack \textrm{Mm}\right\rbrack
% The input-output relations for downsampler in  Z-transform:
% Y\left(z\right)=\frac{1}{M}\sum_{k=0}^{M-1} X\left(z^{1/M} W_M^{-k} \right)
% where W_M =e^{-j2\pi /M}.
% Relation between the discrete-time Fourir transform of the original signal and the discrete-time Fourier transfor of the downsampled signal:
%  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)
%
%% Example 2.1
% This example illustrates the effects of down-sampling in the frequency domain. A signal with triangular shaped spectrum is examined for down-sampling factors M = 2, 3, and 4.

clear all, close all

% Generating the original signal x

F = [0,0.1,0.46,1]; A=[0,1,0,0];  % Setting the input parameters for fir2
x = fir2(256,F,A); % Generating the original signal 'x'
X = fft(x,1024);     % Computing the spectrum of the original signal
f = 0:1/512:(512-1)/512; % Normalized frequencies
figure 
%subplot(4,1,1)
plot(f,abs(X(1:512)))
text(0.7,0.7,'original')
ylabel('|X(e^j^\omega)|'),xlabel('\omega/\pi')
title('Spectrum of the original signal')

% Downsampling-by-2

M = 2;  % Down-sampling factor
y2 = x(1:M:256);   % Down-sampling 
Y2 = fft(y2,1024);   % Computing the spectrum of the down-sampled signal
figure
plot(f,abs(Y2(1:512)))
text(0.6,0.35,'downsampled, M=2')
ylabel('|Y_2(e^j^\omega)|'),xlabel('\omega/\pi')
title('Spectrum of down-sampled signal, M=2')

% Downsampling-by-3

M = 3;  % Down-sampling factor
y3 = x(1:M:256);   % Down-sampling
Y3 = fft(y3,1024);  % Computing the spectrum of the down-sampled signal
figure
plot(f,abs(Y3(1:512)))
text(0.5,0.3,'downsampled, M=3')
ylabel('|Y_3(e^j^\omega)|'),xlabel('\omega/\pi')
title('Effect of aliasing, M=3')

% Downsampling-by-4

M = 4;   % Down-sampling factor
y4=x(1:M:256);   % Down-sampling
Y4 = fft(y4,1024);  % Computing the spectrum of the down-sampled signal
figure
plot(f,abs(Y4(1:512)))
text(0.5,0.2,'downsampled, M=4')
ylabel('|Y_3(e^j^\omega)|'),xlabel('\omega/\pi')
title('Effect of aliasing, M=4')
disp('END OF EXAMPLE 2.1')

%
%% 2.2 Frequency-Domain Characterisation of Up-Sampling
%
% Remainder
%
%
% Factor-of-L upsampler
% The input-output relation of the factor-of-L upsampler:
% y\left\lbrack m\right\rbrack =\left\lbrace \begin{array}{ll}
% x\left\lbrack m/L\right\rbrack , & m=0,\pm L,\pm 2L,\cdots \\
% 0, & \textrm{otherwise}
% \end{array}\right.
% The z-transform of the signal at the output of the factor-of-L upsampler is given by:
% Y\left(z\right)=X\left(z^L \right)
%% Example 2.2
% In this example, we apply the up-sampling with factors L = 2, L = 3,and L = 4 to the original signal x whose spectrum is of the triangular shape.

clear all, close all 

% Generating the input signal

F = [0,0.2,0.9,1]; A = [0,1,0,0];  % Setting the input parameters for fir2
x = fir2(128,F,A); % Generating the original signal ‘x’
X = fft(x,1024);   % Computing the spectrum of the original signal
f = 0:1/512:(512-1)/512;  % Normalized frequencies
figure 
%subplot(4,1,1),
plot(f,abs(X(1:512)))
ylabel('|X(e^j^\omega)|')
text(0.6,0.6,'original')
title('Spectrum of the original signal')

% Upsampling-by-2

L = 2; % Up-sampling factor
y2 = zeros(1,L*length(x)); 
y2([1:L:length(y2)]) = x; % Up-sampled signal, L=2
Y2 = fft(y2,1024); % Computing the spectrum of the up-sampled signal
plot(f,abs(Y2(1:512)))
ylabel('|Y_2(e^j^\omega)|'),xlabel('\omega/\pi') 
text(0.8,0.5,'L=2')
title('Spectrum of the upsampled signal')

% Upsampling-by-3

L = 3; % Up-sampling factor
y3 = zeros(1,L*length(x)); 
y3([1:L:length(y3)]) = x; % Up-sampled signal, L=3
Y3 = fft(y3,1024);  % Computing the spectrum of the up-sampled signal
plot(f,abs(Y3(1:512)))
ylabel('|Y_3(e^j^\omega)|'),xlabel('\omega/\pi')
text(0.75,0.5,'L=3')
title('Spectrum of the upsampled signal')

% Upsampling-by-4

L = 4; % Up-sampling factor
y4 = zeros(1,L*length(x)); 
y4([1:L:length(y4)]) = x; % up-sampled signal, L=4
Y4 = fft(y4,1024);  % Computing the spectrum of the up-sampled signal
plot(f,abs(Y4(1:512)))
ylabel('|(Y_4(e^j^\omega)|'), xlabel('\omega/\pi')
text(0.7,0.5,'L=4')
title('Spectrum of the upsampled signal')
disp('END OF EXAMPLE 2.2')

%
%% 2.3 Spectrum of Decimated Signal
%
% Remainder
%
%
% Block diagram of a decimator: lowpass filter H(z) and factor-of M downsampler
% The role of H(z) is to bandlimit the spectrum of the original signal to \;\pi /M.
%% Example 2.3
% On the example of a triangular-shaped spectrum, we show and compare the effects of decimation with that produced by downsampling.

clear all, close all

% Generating the original signal x
% Setting the input parameters for fir2

F  =[0,0.1,0.46,1]; A = [0,1,0,0]; 
x1 = fir2(256,F,A);
x2 = 0.01*cos(2*pi*0.35*(0:256)); % Generating the signal components

% Original signal

x = x1 + x2;

% Computing the spectrum of the original signal

X = fft(x,1024); % 
f = 0:1/512:(512-1)/512; % Normalized frequencies
figure 
plot(f,abs(X(1:512)))
ylabel('|X(e^j^\omega)|')
text(0.8,0.8,'original')
title('Spectrum of the original signal')

% Down-sampled signal

M = 2;  % Down-sampling factor
y = x(1:M:256);  % Down-sampling
Y = fft(y,1024);   % Computing the spectrum of the down-sampled signal
figure
plot(f,abs(Y(1:512)))
ylabel('|Y(e^j^\omega)|')
text(0.65,0.5,'downsampled M=2')
title('Spectrum of the downsampled signal')

% Decimation by making use of the  the MATLAB function decimate

yd = decimate(x,M);        %  Decimated signal
Yd = fft(yd,1024);         % Computing the spectrum of the decimated signal 
figure
plot(f,abs(Yd(1:512)))
xlabel('\omega/\pi'),ylabel('|Y_d(e^j^\omega)|')
text(0.6,0.3,'decimated, M=2')
title('Spectrum of the interolated signal')
disp('END OF EXAMPLE 2.3')

%
%% 2.4 Spectrum of Interpolated Signal
%
% Remainder
%
%
% Block diagram of an interpolator: A factor-of M upsampler and a lowpass filter H(z).
% The role of H(z) is to bandlimit the spectrum of the upsampled signal to \;\pi /M, and thus eliminate unwanted spectra produced by upsampling.
%% Example 2.4
% On the example of a triangular-shaped spectrum, we show and compare the effects of interpolation with that produced by upsampling.

clear all, close all

% Generation of the original signal x
% Setting the parameters for fir2

F = [0,0.2,0.9,1]; A = [0,1,0,0]; % Setting the input parameters for fir2
x = fir2(128,F,A); % Generating the original signal ‘x’
X = fft(x,1024);  % Computing the spectrum of the original signal
f = 0:1/512:(512-1)/512; % Normalized frequencies
figure 
plot(f,abs(X(1:512)))
ylabel('|X(e^j^\omega)|'),xlabel('\omega/\pi')
text(0.6,0.7,'original')
title('Spectrum of the original signal')
 

% Upsampling-by-4

L = 4;  % Up-sampling factor
xu = zeros(1,L*length(x)); 
xu([1:L:length(xu)]) = x; % Up-sampled signal
Xu = fft(xu,1024);  % Computing the spectrum of the up-sampled signal
figure
plot(f,abs(Xu(1:512)))
ylabel('|X_u(e^j^\omega)|'),xlabel('\omega/\pi')
text(0.62,0.8,'upsampled L=4')
title('Spectrum of the upsampled signal') 

% Interpolation by making use of the  the MATLAB function interp

y = interp(x,L); % Interpolated signal
Y = fft(y,1024); % Computing the spectrum of the up-sampled signal
figure
plot(f,abs(Y(1:512)))
ylabel('|(Y(e^j^\omega)|'), xlabel('\omega/\pi'),
text(0.5,2,'interpolated L=4')
title('Spectrum of the interpolated signal')

% Time-domain presentation

figure 
subplot(3,1,1), stem(54:74,x(55:75))
ylabel('x[n]')
text(68,0.35,'original')
title('Time Domain Prezentation')
subplot(3,1,2)
stem(L*55-1:L*74-1,xu(L*55-1:L*74-1))
ylabel('x_u[m]')
text(270,0.33,'upsampled L=4')
axis([L*55-1,L*74-1,-0.5,0.5])
subplot(3,1,3)
stem(L*55-1:L*74-1,y(L*55-1:L*74-1))
ylabel('y[m]'), xlabel('Time index')
text(270,0.33,'interpolated L=4')
axis([L*55-1,L*74-1,-0.5,0.5])
disp('END OF EXAMPLE 2.4')

%
%% 2.5 The Six Identities
%
%% First identity
%
%% Second identity
%
%% Third Indentity
%
%% Example 2.5
% Illustration of the Third Identity

clear all, close all
n = 0:60;  % Time index
x = cos(2*pi*0.05*n); % Generating the original signal
h = fir1(10,0.5); % Designing the  filter transfer function H(z) 
hu = upsample(h,4); % Transfer function H(z^4)
y1 = filter(hu,1,x); % Filtering
y = downsample(y1,4); % Down-sampling
m = 0:length(y)-1; % Time index
figure  
subplot(3,2,1), stem(n,x), ylabel('x[n]')
title('Figure (a): x[n],y_1[n],y[m]')
subplot(3,2,3), stem(n,y1), ylabel('y_1[n]')
subplot(3,2,5), stem(m,y), ylabel('y[m]')
xlabel('Time index')
y2 = downsample(x,4); % down-sampling
y = filter(h,1,y2); % filtering
subplot(3,2,2), stem(n,x),ylabel('x[n]')
title('Figure (b): x[n],y_2[m],y[m]')
subplot(3,2,4), stem(m,y2), ylabel('y_2[m]')
subplot(3,2,6), stem(m,y),ylabel('y[m]') 
xlabel('Time index')
disp('The third identity verified')
disp('END OF EAMPLE 2.5')

%
%% Fourth Identity
%
%% Fifth Identity
%
%% Sixth Identity
%
%% Example 2.6
% Illustration of the Sixth Identity

clear all, close all
n = 0:15; % Time index
x = cos(2*pi*0.1*n); % Generating the original signal
h = fir1(10,0.5); % Designing the filter transfer function H(z)
hu = upsample(h,2); % Transfer function H(z^2)
y1 = filter(h,1,x); % Filtering
y = upsample(y1,2); % Up-sampling
m = 0:length(y)-1; % Time index
figure  
subplot(3,2,1), stem(n,x), ylabel('x[n]')
title('Figure (a): x[n],y_1[n],y[m]')
subplot(3,2,3), stem(n,y1),ylabel('y_1[n]') 
subplot(3,2,5), stem(m,y), ylabel ('y[m]')
xlabel('Time index')
axis([0,30,-1,1])
y2 = upsample(x,2); % Up-sampling
y = filter(hu,1,y2); % Filtering
subplot(3,2,2), stem(n,x), ylabel('x[n]')
title('Figure (b): x[n],y_2[m],y[m]')
subplot(3,2,4), stem(m,y2), ylabel('y_2[m]')
axis([0,30,-1,1])
subplot(3,2,6), stem(m,y), ylabel('y[m]')
xlabel('Time index')
axis([0,30,-1,1])
disp('The sixth identity verified')
disp('END OF EXAMPLE 2.6')

%
%% 2.6 Cascading Sampling-Rate Alteration Devices
%
% Remainder
%
%
% The two structures shown in the upper Figure are equivalent under the condition that M and L are relativly prime integers.
% We verify the above statement on the example of sampling-rate conversion with the rational factor L/M = 2/3.
%
%% Example 2.7
% Converting the sampling rate by the rational factor L/M=2/3

clear all, close all
n = 0:15;  % Time index
L = 2; M = 3; % Up-sampling and down-sampling factors
x = cos(2*pi*0.1*n); % Generating the original signal
v1 = upsample(x,L);  % Up-sampling
y = downsample(v1,M);  % Down-sampling
r = 0:length(v1)-1;  % Time index
m = 0:length(y)-1;  % Time index
figure 
subplot(3,2,1), stem(n,x), ylabel('x[n]')
title('Figure (a): x[n],v_1[r],y[m]')
subplot(3,2,3), stem(r,v1), ylabel('v_1[r]')
axis([0,30,-1,1])
subplot(3,2,5), stem(m,y), ylabel('y[m]')
xlabel('Time index')
v2 = downsample(x,M);  % Down-sampling
r = 0:length(v2)-1;  % Time index
y = upsample(v2,L); % Up-sampling
m = 0:length(y)-1;  % Time index
subplot(3,2,2), stem(n,x), ylabel('x[n]')
title('Figure (b): x[n],v_2[r],y[m]')
subplot(3,2,4), stem(r,v2), ylabel('v_2[r]')
axis([0,5,-1,1])
subplot(3,2,6), stem(m,y), ylabel('y[m]')
xlabel('Time index')
axis([0,10,-1,1])
disp(['Equivalence satisfied since L and M are relativly prime.'])
disp('END OF EXAMPLE 2.7')

%
%% 2.7 Sampling-Rate Conversion with the Phase Offset
%
% Remainder
% Sampling process can include some phase offset of k samples, i.e., the starting point could be other than zero.Including the phase offset into the down-sampling opperation, we obtain M different signals
%  \left\lbrace y_k \left\lbrack m\right\rbrack \right\rbrace =\left\lbrace x\left\lbrack M\times m+k\right\rbrack \right\rbrace ,\;\;\;\;\;\;\;k=0,1,2,\dots ,M-1
% where for k = 0, the resulting signal \left\lbrace y_0 \left\lbrack m\right\rbrack \right\rbraceis the down-sampled signal without the phase offset. Each value  of k defines a new signal.
% In the next example, we demonstrate the effects of the phase offset to the spectrum of the decimated signal. We show that when the spectrum of the original signal is bandlimitted to \pi /M, there is no aliasing in the signals produced by M-fold decimation. However, when the spectrum of the original signal exceeds this limit, down-sampling opperation produces the aliasing.
%% Example 2.8
% Down-sampling with the phase offset: frequency domain analysis.

clear all, close all

% Generating the original signal x

F = [0,0.1,0.46,1]; A = [0,1,0,0];  % Setting the input parameters for fir2
x = fir2(256,F,A); % Generating the original signal 'x'
X = fft(x,1024);  % Computing the signal spectrum
f = 0:1/512:(512-1)/512; % Normalizad frequencies
figure (1)
subplot(3,2,1), plot(f,abs(X(1:512))), text(0.8,0.7,'(a)')
ylabel('|X(e^{j\omega})|'), axis([0,1,0,1])

% Down-sampling by 2

M = 2; % Down-sampling factor
y0 = downsample(x,M); % down-sampling without phase offset
Y0 = fft(y0,1024); % Computing the signal spectrum
subplot(3,2,3), plot(f,abs(Y0(1:512))), text(0.8,0.7,'(b)')
ylabel('|Y_0(e^{j\omega})|'), axis([0,1,0,1])
 
y1 = downsample(x,M,1); % Down-sampling with the phase offset
Y1 = fft(y1,1024); % Computing the signal spectrum
subplot(3,2,5), plot(f,abs(Y1(1:512))), text(0.8,0.7,'(c)')
ylabel('|Y_1(e^{j\omega})|'), xlabel('\omega/\pi'), axis([0,1,0,1])

% Generating the new original signal

F = [0,0.1,0.7,1]; A = [0,1,0,0]; % Setting the input parameters for fir2
v = fir2(256,F,A);  % Generating the input signal 'v'
V = fft(v,1024);  % Computing the signal spectrum
f = 0:1/512:(512-1)/512;  % Normalized frequencies
figure 
subplot(3,2,2), plot(f,abs(V(1:512))), text(0.8,0.7,'(d)')
ylabel('|V(e^{j\omega})|'), axis([0,1,0,1])
 

% Down-sampling by two

w0=downsample(v,M); % Down-sampling without phase offset
W0=fft(w0,1024); % Computing the signal spectrum
subplot(3,2,4), plot(f,abs(W0(1:512))), text(0.8,0.7,'(e)')
ylabel('|W_0(e^{j\omega})|'), axis([0,1,0,1])
 
w1 = downsample(v,M,1);  % Down-sampling with the phase offset
W1=fft(w1,1024); % Computing the signal spectrum
subplot(3,2,6), plot(f,abs(W1(1:512))), text(0.8,0.7,'(f)')
ylabel('|W_1(e^{j\omega})|'), xlabel(' \omega/\pi'), axis([0,1,0,1])
disp('END OF EXAMPLE 2.8')

% Comment: The plots demonstrate that the magnitude spectrum of the down sampled signal is independent of the phase offset when the frequency band of the original signal is properly bandlimited to \pi /M, see Figs. (a), (b) and (c).  On the contrary,  the right-hand side shows different spectra of the downsampled signals for k = 0 (Fig. (e)) and k=1 (Fig. (f)). Evidently, the ftequency band of the original signal that exceeds \pi /M, Fig. (d) results in different magnitude spectra of the down-sampled  signals.
%
%% 2. 8 Polyphase Decomposition and Reconstruction
%
% This example performs polyphase decomposition of the signal x od the length L = 64 in four polyphase components x0, x1, x2, x3, and reconstruction of the the original signal from the  polyphase components x0, x1, x2, x3.
% Polyphase components x0, x1, x2, x3  are obtained by down-sampling the original signal x with the phase offset.
% In the process of signal reconstruction, up-sampling with the phase offset over x0, x1, x2, x3 is used.
%% Example 2.9

clear all, close all

% Input signal

n = 0:63;
x = zeros(size(n)); x(11:39) = 0.95.^(1:29); % Generating the sequence ‘x’

% Polyphase down-sampling with the phase offset

x0 = downsample(x,4);
x1 = downsample(x,4,1);
x2 = downsample(x,4,2);
x3 = downsample(x,4,3);

% Up-sampling polyphase components with the phase offset

y0 = upsample(x0,4);
y1 = upsample(x1,4,1);
y2 = upsample(x2,4,2);
y3 = upsample(x3,4,3);
figure (1)
subplot (4,1,1), stem(0:length(x0)-1,x0),ylabel('x_0[m]')
title('Polyphase komponents')
subplot (4,1,2), stem(0:length(x1)-1,x1),ylabel('x_1[m]')
subplot (4,1,3), stem(0:length(x2)-1,x2),ylabel('x_2[m]')
subplot (4,1,4), stem(0:length(x3)-1,x3),ylabel('x_3[m]')
xlabel('Time index m')
 

% Signal reconstruction

y = y0 + y1 + y2 + y3; % Adding together up-sampled polyphase components
figure (2)
subplot(2,1,1),stem(n,x), ylabel('x[n]'), axis([0,63,0,1])
title('Original signal')
subplot(2,1,2),stem(n,y), ylabel('y[n]'), axis([0,63,0,1])
title('Reconstructed signal')
xlabel('Time index n')
disp('END OF EXAMPLE 2.9')

%
%% 2.9 Multistage Systems
%
% Equivalent frequency response of multistage decimator.
% When decimation factor M can be factored into K integers M=M_1 \times M_2 \times \dots M_K, instead of using a single filter and factor-of-M down-sampler the overall decimator can be implemented as a cascade of K decimators as shown bellow.
%
% The transfer function H\left(z\right)of the equivalent single-stage decimation filter can be obtained by applying the third identity to the above implementation scheme:
% H\left(z\right)=H_1 \left(z\right)\timesH_2 \left(z^{M_1 } \right)\times H_3 \left(z^{{M_1 M}_2 } \right)\times \cdots H_K \left(z^{{M_1 M}_2 \cdots M_{K-1} } \right)
%% Example 2.10
% This program computes the equivalent frequency response of the example three-stage factor-of-8 decimator.

clear all, close all
h1 = fir1(30,0.5); % FIR filter H1(z)
[H1,f] = freqz(h1,1,512,2);
figure 
subplot(4,1,1), plot(f,abs(H1)), ylabel('|H_1(e^{j\omega})|')
axis([0,1,0,1.2])
 
M = 2;
h2 = zeros(1,M*length(h1)); 
h2([1:M:length(h2)]) = h1; % H2(z) = H1(z^2)
[H2,f] = freqz(h2,1,512,2);
subplot(4,1,2), plot(f,abs(H2)), ylabel('|H_1(e^{j2\omega})|')
axis([0,1,0,1.2])
 
M = 4;
h3 = zeros(1,M*length(h1)); 
h3([1:M:length(h3)]) = h1; %  H3(z) = H1(z^4)
[H3,f] = freqz(h3,1,512,2);
subplot(4,1,3), plot(f,abs(H3)), ylabel('|H_1(e^{j4\omega})|')
axis([0,1,0,1.2])
 
H = H1.*H2.*H3; % Equivalent filter
subplot(4,1,4), plot(f,abs(H)),ylabel('|(H(e^{j\omega})|') 
xlabel('\omega/\pi'), axis([0,1,0,1.2])
disp('END OF EXAMPLE 2.10')
disp('                    END OF CHAPTER II')

%

