%%
% by Ljiljana Milic
% Supplemental material for Chapter I:
%% 1. Single-Rate Signals and Systems: Background Review
% Table of Contents
% 1. Single-Rate Signals and Systems: Background Review
%      1.1 Discrete Time Fourier Transform (DTFT)
%         Example 1.1
%     1.2 Discrete Fourier Transform (DFT)
%         Example 1.2
%     1.3 Linear Time-Invariant (LTI) System
%         Example 1.3
%%  1.1 Discrete Time Fourier Transform (DTFT)
%
% Remainder
% Definition of the discrete-time Fourier transform (DTFT):
% X\left(e^{j\omega } \right)=\sum_{n=-\infty }^{\infty } x\left\lbrack n\right\rbrack e^{-j\omega n}
%  DTFT for the test signal \left\lbrace x\left\lbrack n\right\rbrack \right\rbraceof a finite length L
% X\left(e^{j\omega } \right)=\sum_{n=0}^L x\left\lbrack n\right\rbrack e^{-j\omega n}
%% Example 1.1
% In this example, the spectrum of the triangular signal \left\lbrace x\left\lbrack n\right\rbrack \right\rbrace is computed using function freqz from the Signal Processing Toolbox.

clear all, close all

% Introducing the test signal \left\lbrace x\left\lbrack n\right\rbrack \right\rbrace

L = 13;
x = triang(13);
L = length(x); % length of the sequence
N = 256;  % number of DTFT points        

% Computing Discrete Time Fourier Transform X\left(e^{j\omega } \right) using Matlab function freqz

[X,w] = freqz(x,1,N); % computation of the signal spectrum 
mag = abs(X); % magnitude spectrum
phase = angle(X); % phase spectrum

% Displaying the results =

figure
stem(0:L-1,x), 
title('Signal {x[n]}')
xlabel('Time index n'), ylabel('x[n]') 
figure
plot(w/pi,mag), 
title('Magnitude Spectrum')
xlabel('Normalized frequency \omega/\pi'), ylabel('|X(e^{j\omega})|')
figure
plot(w/pi,unwrap(phase)),
title('Phase Spectrum')
xlabel('Normalized frequency \omega/\pi'), ylabel('\phi(\omega)')
disp('END OF EXAMPLE 1.1')

%
%% 1.2 Discrete Fourier Transform (DFT)
%
% Remainder
% Definition of the discrete Fourier transform (DFT):
% X\left\lbrack k\right\rbrack ={X\left(e^{|j\omega } \right)}_{|\omega =2\pi k/N} =\sum_{n=0}^{L-1} x\left\lbrack n\right\rbrack e^{-j\pi \textrm{kn}/N},   0\le k\le N-1
% where N\ge L
% With substitution W_N =e^{-\textrm{j2}\pi /N} , DFT is presented in the form:
%  X\left\lbrack k\right\rbrack =\sum_{n=0}^{L-1} x\left\lbrack n\right\rbrack W_N^{-\textrm{nk}},   0\le k\le N-1
%% Example 1.2
% In this example, we use function fft to compute DFT sequence \left\lbrace X\left\lbrack k\right\rbrack \right\rbrace for the test signal \left\lbrace x\left\lbrack n\right\rbrack \right\rbracecomposed of three sinusoidal components .

clear all close all
n = 0:63; % Time index n

% Introducing test signal x[n]

x=sin(2*pi*4*n/64)+0.6*sin(2*pi*8*n/64) + 0.8*sin(2*pi*18.5*n/64);
figure 
stem(n,x)
xlabel('Time index n'), ylabel('x[n]')
title('Signal x[n]')
axis([0,63,-3,3])

% Computing Discrete Fourier Trnsform X[k]

X = fft(x); % Computation of DFT
k = n;     % Frequency index k
figure
stem(k,abs(X))
xlabel('Frequency index k'), ylabel('|(X[k])|')
title('DFT sequence')
axis([0,63,0,40])
disp('END OF EXAMPLE 1.2')

%
%% 1.3 Linear Time-Invariant (LTI) System
% Computation of frequency response and pole-zero plot
% Remainder
%
% Transfer function of an LTI system H\left(z\right)
% H\left(z\right)=\frac{\sum_{k=0}^M b_k z^{-k} }{1+\sum_{k=1}^N a_k z^{-k} }
% Frequency response
% H\left(e^{j\omega } \right)=H{\left(z\right)}_{|z=e^{j\omega } }
% Magnitude response: |H\left(e^{j\omega } \right)|
% Phase response \phi \left(\omega \right)=\arg \left\lbrack H\left(e^{j\omega } \right)\right\rbrack \;
% Group delay \tau \left(\omega \right)=-\frac{d\phi \left(\omega \right)}{d\omega }
%% Example 1.3
% LTI system example: Design and analysis of the 5th order Chebyshev filter by using function cheby1 for design.

clear all, close all

% Designing the transfer function H\left(z\right)

[B,A] = cheby1(5,1,0.4); % 5th order Chebyshev filter

% Computing frequency response H\left(e^{j\omega } \right)

[H,f] = freqz(B,A,250,2); % Frequency response
Mag=abs(H); % Magnitude response
Phase = unwrap(angle(H));  % Phase response
[Gd,f] = grpdelay(B,A,250,2); % Group delay 

% Displaying the results

figure 
zplane(B,A),text(1,0.9,'(a)') % Pole-zaro plot
title('Pole-zero locations')
figure
plot(f,Mag), axis([0,1,0,1.1]),text(0.8,0.97,'(b)')
xlabel('\omega/\pi'), ylabel('Magnitude')
title('Magnitude response')
plot(f,Phase), axis([0,1,-8,0]),text(0.8,-1,'(c)')
xlabel(' \omega/\pi'), ylabel('Phase, rad')
title('Phase response')
plot(f,Gd), axis([0,1,0,15])
xlabel('\omega/\pi'), ylabel('Group delay, samples'),text(0.8,13,'(d)')
title('Group delay')
disp('END OF EXAMPLE 1.3')
disp('                     END OF CHAPTER I')

%

