function g = czt(x, k, w, a)
%CZT	Chirp-z transform.
%	G = CZT(X, K, W, A) returns:
%	G   the K-element chirp-z transform of data X,
%	X   the input matrix of signal columns
%	K   the length of the chirp z-transform
%	W   is the ratio W between points along the
%	complex-plane spiral contour of interest
%	A   the complex starting point on that contour

% References:
%  Oppenheim, A.V. & R.W. Schafer, Discrete-Time Signal
%   Processing,  Prentice-Hall, pp. 623-628, 1989.
%  Rabiner, L.R. and B. Gold, Theory and Application of
%   Digital Signal Processing, Prentice-Hall, Englewood
%   Cliffs, New Jersey, pp. 393-399, 1975.

% Charles R. Denham, MathWorks, 1990.
% modified: jMcClellan, 7/25/90.
% modified: cDenham, 8/15/90.

[m, n] = size(x); oldm = m;
if m == 1, x = x(:); [m, n] = size(x); end

if nargin < 2, k = length(x); end
if nargin < 3, w = exp(-sqrt(-1) .* 2 .* pi ./ k); end
if nargin < 4, a = 1; end

%------- Length for power-of-two fft.

nfft = 1;
while nfft < (m + k - 1), nfft = 2 .* nfft; end

%------- Premultiply data.

kk = ( (-m+1):(k-1) ).';
kk2 = (kk .^ 2) ./ 2;
ww = w .^ (kk2);   % <----- Chirp filter is 1./ww
aa = a .^ (  -kk( m:(2*m-1) )   );
y = x .* [   ( aa.*ww( m:(2*m-1) )  ) * ones(1, n)   ];

%------- Fast convolution via FFT.

fy = fft(  y, nfft );
fv = fft( 1 ./ ww, nfft );   % <----- Chirp filter.
fy = fy .* ( fv * ones(1, n) );
g  = ifft( fy );

%------- Final multiply.

g = g( m:(m+k-1), : ) .* [ ww( m:(m+k-1) )  *  ones(1, n)   ];

if oldm == 1, g = g.'; end
