home *** CD-ROM | disk | FTP | other *** search
/ Liren Large Software Subsidy 10 / 10.iso / l / l455 / 8.ddi / SIGNAL.DI$ / CZT.M < prev    next >
Encoding:
Text File  |  1993-03-11  |  1.5 KB  |  54 lines
  1. function g = czt(x, k, w, a)
  2. %CZT    Chirp-z transform.
  3. %    G = CZT(X, K, W, A) returns:
  4. %    G   the K-element chirp-z transform of data X,
  5. %    X   the input matrix of signal columns
  6. %    K   the length of the chirp z-transform
  7. %    W   is the ratio W between points along the
  8. %    complex-plane spiral contour of interest
  9. %    A   the complex starting point on that contour
  10.  
  11. % References:
  12. %  Oppenheim, A.V. & R.W. Schafer, Discrete-Time Signal
  13. %   Processing,  Prentice-Hall, pp. 623-628, 1989.
  14. %  Rabiner, L.R. and B. Gold, Theory and Application of
  15. %   Digital Signal Processing, Prentice-Hall, Englewood
  16. %   Cliffs, New Jersey, pp. 393-399, 1975.
  17.  
  18. % Charles R. Denham, MathWorks, 1990.
  19. % modified: jMcClellan, 7/25/90.
  20. % modified: cDenham, 8/15/90.
  21.  
  22. [m, n] = size(x); oldm = m;
  23. if m == 1, x = x(:); [m, n] = size(x); end
  24.  
  25. if nargin < 2, k = length(x); end
  26. if nargin < 3, w = exp(-sqrt(-1) .* 2 .* pi ./ k); end
  27. if nargin < 4, a = 1; end
  28.  
  29. %------- Length for power-of-two fft.
  30.  
  31. nfft = 1;
  32. while nfft < (m + k - 1), nfft = 2 .* nfft; end
  33.  
  34. %------- Premultiply data.
  35.  
  36. kk = ( (-m+1):(k-1) ).';
  37. kk2 = (kk .^ 2) ./ 2;
  38. ww = w .^ (kk2);   % <----- Chirp filter is 1./ww
  39. aa = a .^ (  -kk( m:(2*m-1) )   );
  40. y = x .* [   ( aa.*ww( m:(2*m-1) )  ) * ones(1, n)   ];
  41.  
  42. %------- Fast convolution via FFT.
  43.  
  44. fy = fft(  y, nfft );
  45. fv = fft( 1 ./ ww, nfft );   % <----- Chirp filter.
  46. fy = fy .* ( fv * ones(1, n) );
  47. g  = ifft( fy );
  48.  
  49. %------- Final multiply.
  50.  
  51. g = g( m:(m+k-1), : ) .* [ ww( m:(m+k-1) )  *  ones(1, n)   ];
  52.  
  53. if oldm == 1, g = g.'; end
  54.