%========================================================
% Amine Faridi (29061780)
% Project # 5
% EE416
%========================================================
%========================================================
% a)-
%========================================================

% 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 .
% Strictly speaking, the `chirp transform' finds these samples along 
% the unit circle, while the `chirp z-transform' algorithm (CTA) refers to samples 
% along a general contour (which, of course, may be the unit circle).
% From the discussion in the class:
% If we have  x[n] as a N-point sequence and, 
%  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)   when  -(N-1) <= n <= M-1
%                      = 0                     otherwise

%========================================================
% b)-
%========================================================

% The advantage of the chirp z-transform is that it enables 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. In other words, if the 
% formant frequency of interest is in the range 200 Hz to 900 Hz, then samples of 
% the z-transform may be calculated specifically in this region. In contrast, the 
% frequency range of the DFT is strictly related to the sampling frequency. 
% Furthermore, the chirp z-transform allows the calculations of arbitrary number of 
% samples along this contour, which arbitrarily reduces the error in frequency 
% representation. 
% As discussed previously for the DFT case, if a frequency representation of higher 
% resolution is desired, the only recourse is to zero-pad the signal, as previously 
% done in project 3. 
% Finally, the chirp z-transform allows one to leave the unit circle and calculate 
% samples in the z-plane, 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.




%========================================================
% c)-
%========================================================

fs=1000;		% The stop frequency
f=15;			% Frequency 15Hz
w=2*pi*f;		% To find angular frequency
dw=2*pi*3/4;	% To increment frequency
t=0:1/fs:1.5;	% The time range



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

h=fir1(50,[0.015 0.035],kaiser(51,0.5));	% The FIR1 filter using
                                            % the Kaiser window
y=conv(x,h);								% The filter output
Yz=freqz(y);								% The Z-transform of 
                                            % the filter output

figure;													
plot(abs(Yz));								% Plot Z-transform Yz
title('Plot of Z-Transform using normal FT function');		
xlabel('Freq');										
ylabel('amplitude');										

f1=7.5;
f2=17.5;

m=256;									% The number of elements									
w=exp(-j*2*pi*(f2-f1)/(m*fs));			% The complex ratio
a=exp(j*2*pi*f1/fs);					% The complex starting point

z=czt(y,m,w,a);						    % CZT function
fz=((0:length(z)-1)'*(f2-f1)/length(z))+f1;% frequency range

figure;												
plot(fz,abs(z));						% Ploting Chirp Z-transform 
axis([f1 f2 0 4000]);						

clear all;

% For other frequencies, and parameters

fs=3000;		
f=20;			
w=2*pi*f;		
dw=2*pi*3/4;	
t=0:1/fs:1.5;	



x=sin(w*t)+2*sin((w+dw)*t)+3*cos((w-dw)*t);	
Xz=freqz(x);								

h=fir1(50,[0.015 0.035],kaiser(51,0.5));	                                        
y=conv(x,h);								
Yz=freqz(y);								
                                           

figure;													
plot(abs(Yz));								
title('Plot of Z-Transform using normal FT function');		
xlabel('Freq');										
ylabel('amplitude');										

f1=10.5;
f2=27.5;

m=256;																		
w=exp(-j*2*pi*(f2-f1)/(m*fs));			
a=exp(j*2*pi*f1/fs);						

z=czt(y,m,w,a);						
fz=((0:length(z)-1)'*(f2-f1)/length(z))+f1;

figure;												
plot(fz,abs(z));							 
axis([f1 f2 0 8000]);						

% For another frequencies and parameters:

clear all;

fs=500;		
f=10;			
w=2*pi*f;		
dw=2*pi*3/4;	
t=0:1/fs:1.5;	



x=sin(w*t)+2*sin((w+dw)*t)+3*cos((w-dw)*t);	
Xz=freqz(x);								

h=fir1(70,[0.025 0.065],kaiser(71,0.5));	
                                           
y=conv(x,h);								
Yz=freqz(y);								 
                                            

figure;													
plot(abs(Yz));								
title('Plot of Z-Transform using normal FT function');		
xlabel('Freq');										
ylabel('amplitude');										

f1=2.5;
f2=17.5;

m=512;																		
w=exp(-j*2*pi*(f2-f1)/(m*fs));			
a=exp(j*2*pi*f1/fs);						

z=czt(y,m,w,a);						
fz=((0:length(z)-1)'*(f2-f1)/length(z))+f1;

figure;												
plot(fz,abs(z));							
axis([f1 f2 0 4000]);						
