%Vedrana Sarajlic
%ID#2730-3757
%EE416 Project 5
%12/3/01

%The Chirp Transform Algorithm (CTA) is another algorithm based on expressing the DFT
%as a convolution; however, it is not optimal in minimizing arithmetic complexity.
%It can be used to compute any set of equally spaced samples of the Fourier transform
%on the unit circle. 
%The Chirp Z-Transform (CZT) is a more general form of the CTA; it is a way by which one
%can determine samples of the z-transform equally spaced along a contour in the z-plane.

%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 contrast, the frequency range of the DFT
%is strictly related to the sampling frequency.  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.  In addition, the CZT allows one to leave
%the unit circle and calculate samples in the z-plane, while the DFT is limited to the 
%unit circle.

clear all;

%let
t=0:1/1000:2;
f=15;
w=2*pi*f;
dw=2*pi*0.5;

x=sin(w*t)+2*sin((w+dw)*t)+3*cos((w-dw)*t);
figure(1)
plot(x);
title('Signal x=sin(wt)+2sin((w+dw)t)+3cos((w-dw)t)');

%convolve the signal with an FIR filter
h=fir1(32,.5,kaiser(33,4));
y=conv(x,h);
figure(2);

freqz(y);
title('Frequency response');

%Chirp-Z transform of the signal

M=length(x); %parameters (default)
W=exp(-j*2*pi/M);
A=1;
fs=1000;
f1=10;
f2=1050;

ch=czt(x,M,W,A);
figure(3);
fy=(0:length(y)-1)'*1000/length(y);
fch=((0:length(ch)-1)'*(f2-f1)/length(ch))+f1;

subplot(211);
plot(fy(1:500),abs(fft(y(1:500))));axis([1 250 0 50])
title('FFT');
subplot(212);
plot(fch,abs(ch));axis([f1 f2 0 20]);
title('CZT');



%change the values of the parameters
fs=1000;
f1=100;
f2=1000; %in hertz
M1=1024;
W1=exp(-j*2*pi*(f2-f1)/(M1*fs));
A1=exp(j*2*pi*f1/fs);
y1=fft(y,1000);
z=czt(x,M1,W1,A1);
fy1=(0:length(y1)-1)'*1000/length(y1);
fz=((0:length(z)-1)'*(f2-f1)/length(z))+f1;
figure(4);
subplot(211);
plot(fy1(1:500),abs(y1(1:500)));axis([1 300 0 20])
title('FFT');
subplot(212);
plot(fz,abs(z));axis([f1 f2 0 10]);
title('CZT');


figure(5);
M2=500;
W2=exp(-j*2*pi*(fs-f1)/(M2*fs));
A2=exp(j*2*pi*f1/fs);
y2=fft(y);
ch1=czt(x,M2,W2,A2);
fy2=(0:length(y2)-1)'*1000/length(y2);
fch1=((0:length(ch1)-1)'*(f2-f1)/length(ch1))+f1;
subplot(211);
plot(fy2(1:500),abs(y2(1:500)));axis([1 250 0 30])
title('FFT');
subplot(212);
plot(fch1,abs(ch1));axis([100 500 0 10]);
title('CZT');

