
%  **************************************************************************

%  								MATLAB    PROJECT    5

%  This file will analyze the concept of Chirp Z-transform and application

%  **************************************************************************


%  ------------------------------------------------------------------------
%  (1) Explain Chirp Z-transform (CZT) and Chirp Transform Algorithm (CTA)
%  ------------------------------------------------------------------------

%  Chirp Transform Algorithm (CTA) is another expression of DFT as a convolution. This
%  algorithm does not reduce arithmetic complexity, but well suited for applications
%  that involved coonvolution of a fixed, prespecified impulse response. CTA does not
%  go around the extire unit circle, unlike FFT, so it can compute any set of equally
%  spaced samples of Fourier Transform.

%  The M samples are spaced equally at frequencies
%                    wk = w0 + k(Dw),          k = 0,1...,M-1

%  The Chirp Transform Algorithm (CTA) is :

%           X(ejwn) = (sum g[k]W^((-(n-k)^2)/2)) W^(n^2/2)

%  where g[n] = x[n]e(-jwn)w^(n^2/2)

%  Thus, the impulse response is defined as

%                 h(n) = w^(-n^2/2)       for  -(N-1) <= n <= M-1
%                      = 0                     otherwise


%  ------------------------------------------------------------------------
%  (2) Enumerate key advantages and applications of Chirp Z-Transform (CZT)
%  ------------------------------------------------------------------------


%  2)	The advantage of the chirp z-transform is that it enables 
%     (a)	one to calculate the samples of the z-transform equally spaced over an arc or a 
%           spiral contour with an arbitrary starting point and arbitrary frequency range; 
%     (b)	In contrast, the frequency range of the DFT is strictly related to the sampling frequency. 
%     (c)	Furthermore, the chirp z-transform allows one to calculate an arbitrary number of samples 
%           along this contour, which arbitrarily reduces the error in frequency representation. 
%     (d)	As discussed previously for the DFT case, if one requires a frequency representation 
%           of higher resolution, the only recourse is to zero-pad the signal (i.e., the number of 
%           sample points in the DFT is the same as the number of points in the signal).
%     (e)	Lastly, the chirp z-transform allows one to leave the unit circle and calculate
%           samples in the  , while the DFT is limited to the unit circle. This is particularly 
%           useful for cases in which the frequencies corresponding to poles inside the unit circle
%           cannot be distinguished along the unit circle contour. 

%     However, there is a disadvantage of the chirp-z transform is that it may be slower than 
%     the FFT algorithm, depending on the number of sample points


%  -----------------------------------------
%  (3) Use Matlab to show advatage of CZT
%  -----------------------------------------

fs=1000;														% Set stop frequency
f=15;															% Set frequency to be 60Hz
w=2*pi*f;													% Find angular frequency
dw=2*pi*3/4;												% Set frequency increment
t=0:1/fs:1.5;												% Set time ranges

x=sin(w*t)+2*sin((w+dw)*t)+3*cos((w-dw)*t);		% Define signal
Xz=freqz(x);												% Find Z-Transform of x

h1=fir1(50,[0.015 0.035],kaiser(51,0.5));			% Define filter
y1=conv(x,h1);												% Find filter output
Yz=freqz(y1);												% Find Z-transform of y1
	
figure(1);													% Set plot 1
plot(abs(Yz));												% Plot Z-transform Yz
title('Plot of Z-Transform using normal FT function');		% Label title
xlabel('frequency');										% Label x-axis
ylabel('amplitude');										% Label y-axis

f1=7.5;
f2=17.5;
m1=512;														% Set number of elements									
w1=exp(-j*2*pi*(f2-f1)/(m1*fs));						% Set complex ratio
a1=exp(j*2*pi*f1/fs);									% Set complex starting point

z1=czt(y1,m1,w1,a1);										% Compute using CZT function
fz1=((0:length(z1)-1)'*(f2-f1)/length(z1))+f1;	% Set frequency range

figure(2);													% Set plot 2
plot(fz1,abs(z1));										% Plot the Chirp Z-transform computed
axis([f1 f2 0 4000]);									% Set axis size
title('Plot of Chirp Z-Transform for fs=1000Hz');			% Label title
xlabel('amplitude');										% Label x-axis
ylabel('frequency');										% Label y-axis


fs=2000;														% Set stop frequency
ft2=20;														% Set frequency to be 20Hz
w2=2*pi*ft2;												% Find angular frequency
dw=2*pi*3/4;												% Set frequency increment
t=0:1/fs:1.5;												% Set time ranges

x2=sin(w2*t)+2*sin((w2+dw)*t)+3*cos((w2-dw)*t);	% Define signal
Xz=freqz(x2);												% Find Z-Transform of x

h2=fir1(60,[0.035 0.055],kaiser(61,0.6));			% Define filter
y2=conv(x2,h2);											% Find filter output
Yz=freqz(y2);												% Find Z-transform of y2
	

f2=17.5;
f3=27.5;
m2=1028;														% Set number of elements									
w2=exp(-j*2*pi*(f3-f2)/(m2*fs));						% Set complex ratio
a2=exp(j*2*pi*f2/fs);									% Set complex starting point

z2=czt(y2,m2,w2,a2);										% Compute using CZT function
fz2=((0:length(z2)-1)'*(f3-f2)/length(z2))+f2;	% Set frequency range

figure(3);													% Set plot 2
plot(fz2,abs(z2));										% Plot the Chirp Z-transform computed
axis([f2 f3 0 2500]);									% Set axis size
title('Plot of Chirp Z-Transform for fs=2000Hz');				% Label title
xlabel('amplitude');										% Label x-axis
ylabel('frequency');										% Label y-axis



fs=3500;														% Set stop frequency
ft3=25;														% Set frequency to be 25Hz
w3=2*pi*ft3;												% Find angular frequency
dw=2*pi*3/4;												% Set frequency increment
t=0:1/fs:1.5;												% Set time ranges

x3=sin(w3*t)+2*sin((w3+dw)*t)+3*cos((w3-dw)*t);	% Define signal
Xz=freqz(x3);												% Find Z-Transform of x

h3=fir1(70,[0.055 0.075],kaiser(71,0.7));			% Define filter
y3=conv(x3,h3);											% Find filter output
Yz=freqz(y3);												% Find Z-transform of y1
	
f3=19;
f4=29;
m3=1028;														% Set number of elements									
w3=exp(-j*2*pi*(f4-f3)/(m3*fs));						% Set complex ratio
a3=exp(j*2*pi*f3/fs);									% Set complex starting point

z3=czt(y3,m3,w3,a3);										% Compute using CZT function
fz3=((0:length(z3)-1)'*(f4-f3)/length(z3))+f3;	% Set frequency range

figure(4);													% Set plot 2
plot(fz3,abs(z3));										% Plot the Chirp Z-transform computed
axis([f3 f4 0 2000]);									% Set axis size
title('Plot of Chirp Z-Transform for fs=3500Hz');				% Label title
xlabel('amplitude');										% Label x-axis
ylabel('frequency');										% Label y-axis



%  The plots show that the plots of Chirp Z-Transform are have finite sequence
%  length, because it does not goes arount the unit circle like DFT. CZT plot
%  will approach DFT when the change becomes a period with resprct to N samples 
%  (or 2pi/N ). 


% ----------------------------------------------

% NAME    : FOONG CHAN TOONG

% UB NO.  : 2764 - 1508

% CLASS   : EE 416

% PROJECT : MATLAB PROJECT 5

% ----------------------------------------------
