c10ex16.m

为了更好的便于增加交流

%	Implementation of the chirp-z transform of a pulse signal and
%	comparison with the corresponding results for the FFT - Example 10-16
%
clf
theta_0 = input('Enter starting point, theta_0, on unit circle in radians ');
phi_0 = input('Enter frequency resolution, phi_0, around unit circle in radians ');
N = input('Enter total number of points for pulse signal ');
N_pulse = input('Enter number of nonzero points for pulse signal ');
M = input('Enter number of frequency cells for chirp-z transform ');
W = exp(-j*phi_0);
A = exp(j*theta_0);
n = 0:N-1;
m = 0:M-1;
x = zeros(size(1:N));
for nn = 1:N_pulse
	x(nn) = 1;
end
k = 1:M;
T = A.*W.^(-(k-1));
Y = ones(size(T));	% Form matrix of (AW^(-1))^(k*n); k = columns & n = rows
for nn = 2:N
	Y = [Y; T.^(nn-1)];
end
X_chirp = x*Y;		% Use matrix multiply to compute chirp-z sum
X_FFT = fft(x);		% Compute FFT for comparison
%	Shift chirp-z and set up axes in plot to make direct comparison
%	between chirp-z and FFT results.
NN = 2*pi/phi_0-1;
zz = zeros(size(1:theta_0*(NN+1)/(2*pi)));
X_chirp = [zz X_chirp];
L = length(X_chirp);
l = 0:L-1;
subplot(3,1,1),stem(n,x),axis([0 N-1 0 1]),xlabel('n'),ylabel('x(n)')
subplot(3,1,2),stem(l,abs(X_chirp)),axis([0 NN-1 0 N_pulse]),xlabel('k_chirp'),...
	ylabel('Magn. of chirp-z'),...
	text(NN/2,.7*N_pulse,['Resolution = ',num2str(phi_0),' radians'])
subplot(3,1,3), stem(n,abs(X_FFT)), axis([0 N-1 0 N_pulse]),xlabel('k_FFT'),...
	ylabel('Magn. of FFT')