 % ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
%  EE 416
%  Hok-Lim Li (29045906)
%  Matlab Project #5 %  "Chirp Z-transform and application" 
%  **************************************************************************   
%  (1) Explain Chirp Z-transform (CZT) and Chirp Transform Algorithm (CTA) 
%  1)	The chirp z-transform algorithm is a means by which one can determine samples of 
%  the z-transform equally spaced along a contour in the z-plane. The 'chirp transform' 
%  finds these samples along the unit circle, while the `chirp z-transform' algorithm 
%  refers to samples along a general contour.  
%  Also, 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 and prespecified impulse response. CTA can compute any set
%  of equally spaced samples of Fourier Transform. 
%  The M samples are spaced equally at frequencies wk = w0 + k(Dw) when 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). So, the impulse response          %          {       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 advantages of the chirp z-transform are:  %     (a)	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)	The frequency range of the DFT is strictly related to the sampling frequency.  %     (c)	The CZT allows one to calculate an arbitrary number of samples  %           along this contour, which arbitrarily reduces the error in frequency representation.  %     (d)	CZT provides a frequency representation of higher resolution. %     (e)	The CZT 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.   %     There maybe 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  stopf=1000;													% Set stop frequency 

ft1=20;														% Set frequency to be 20Hz 
w=2*pi*ft1;													% Define angular frequency 
dw=2*pi*3/4;												% Set frequency increment 
t=0:1/stopf:1.5;											% Set time ranges  
X1=sin(w*t)+2*sin((w+dw)*t)+3*cos((w-dw)*t);		% Define signal 
Xz=freqz(X1);												% Find Z-Transform of x  
H1=fir1(50,[0.015 0.035],kaiser(51,0.5));			% Define filter 
y1=conv(X1,H1);											% Find filter output 
Yz=freqz(y1);												% Find Z-transform of y1 	 
figure(1);	
   
plot(abs(Yz));												% Plot Z-transform Yz 
title('Z-Transform with normal FT function');	 
xlabel('f');							

ylabel('amplitude');										 

F1=7.5; 
F2=17.5; 
m1=512;														% Set number of elements									 
w1=exp(-j*2*pi*(F2-F1)/(m1*stopf));					% Set complex ratio 
a1=exp(j*2*pi*F1/stopf);								% 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('Chirp Z-Transform for 
=100Hz');		

xlabel('amplitude');										

ylabel('f');												  

stopf=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/stopf: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*stopf));					% Set complex ratio 
a2=exp(j*2*pi*F2/stopf);								% 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('Chirp Z-Transform stopf=2000Hz');			

xlabel('amplitude');										 
ylabel('f');										 
stopf=3500;													% Set stop frequency 
ft3=25;														% Set frequency to be 25Hz 
w3=2*pi*ft3;												% Angular frequency 
dw=2*pi*3/4;												% Set frequency increment 
t=0:1/stopf: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 element numbers									 
w3=exp(-j*2*pi*(F4-F3)/(m3*stopf));					% Set complex ratio 
a3=exp(j*2*pi*F3/stopf);								% Define 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);													

plot(fz3,abs(z3));										% Plot the Chirp Z-transform computed 
axis([F3 F4 0 1600]);									% Set axis size 
title('CZT for stopf=3500Hz');						

xlabel('Amplitude');										

ylabel('f');												 