[OctDev] (no subject)

[<dg...@um...>]

(It seems the octave-dev list is very quiet, so I'm also sending this to the
octave-help list.)

The czt function was not working, so I implemented it correctly.  At a
minimum, a correct czt should pass this test:

x=3Drand(1000,1);
y1=3Dfft(x);
y2=3Dczt(x);
max(abs(y1-y2)) # should be a very small number

Also fixed the documentation.

Will someone put this in the cvs?  Thanks.

- DG

## Copyright (C) 2004 Daniel Gunyan
##
## This program is free software; you can redistribute it and/or modify
## it under the terms of the GNU General Public License as published by
## the Free Software Foundation; either version 2 of the License, or
## (at your option) any later version.
##
## This program is distributed in the hope that it will be useful,
## but WITHOUT ANY WARRANTY; without even the implied warranty of
## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
## GNU General Public License for more details.
##
## You should have received a copy of the GNU General Public License
## along with this program; if not, write to the Free Software
## Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA

## usage y=3Dczt(x, m, w, a)
##
## Chirp z-transform.  Compute the frequency response starting at a and
## stepping by w for m steps.  a is a point in the complex plane, and
## w is the ratio between points in each step (i.e., radius increases
## exponentially, and angle increases linearly).
##
## To evaluate the frequency response for the range f1 to f2 in a signal
## with sampling frequency Fs, use the following:
##     m =3D 32;                               ## number of points desired
##     w =3D exp(-j*2*pi*(f2-f1)/((m-1)*Fs));  ## freq. step of f2-f1/m
##     a =3D exp(j*2*pi*f1/Fs);                ## starting at frequency f1
##     y =3D czt(x, m, w, a);
##
## If you don't specify them, then the parameters default to a fourier=20
## transform:
##     m=3Dlength(x), w=3Dexp(-j*2*pi/m), a=3D1
##
## If x is a matrix, the transform will be performed column-by-column.

## Algorithm (based on Oppenheim and Schafer, "Discrete-Time Signal
## Processing", pp. 623-628):
##   make chirp of length -N+1 to max(N-1,M-1)
##     chirp =3D> w^([-N+1:max(N-1,M-1)]^2/2)
##   multiply x by chirped a and by N-elements of chirp, and call it g
##   convolve g with inverse chirp, and call it gg
##     pad ffts so that multiplication works
##     ifft(fft(g)*fft(1/chirp))
##   multiply gg by M-elements of chirp and call it done

function y =3D czt(x, m, w, a)
  if nargin < 1 || nargin > 4, usage("y=3Dczt(x, m, w, a)"); endif

  [row, col] =3D size(x);
  if row =3D=3D 1, x =3D x(:); col =3D 1; endif

  if nargin < 2 || isempty(m), m =3D length(x(:,1)); endif
  if length(m) > 1, error("czt: m must be a single element\n"); endif
  if nargin < 3 || isempty(w), w =3D exp(-2*j*pi/m); endif
  if nargin < 4 || isempty(a), a =3D 1; endif
  if length(w) > 1, error("czt: w must be a single element\n"); endif
  if length(a) > 1, error("czt: a must be a single element\n"); endif

  ## indexing to make the statements a little more compact
  n =3D length(x(:,1));
  N =3D [0:n-1]'+n;
  NM =3D [-(n-1):(m-1)]'+n;
  M =3D [0:m-1]'+n;

  nfft =3D 2^nextpow2(n+m-1); # fft pad
  W2 =3D w.^(([-(n-1):max(m-1,n-1)]'.^2)/2); # chirp

  for idx =3D 1:col
    fg =3D fft(x(:,idx).*(a.^-(N-n)).*W2(N), nfft);
    fw =3D fft(1./W2(NM), nfft);
    gg =3D ifft(fg.*fw, nfft);

    y(:,idx) =3D gg(M).*W2(M);
  endfor

  if row =3D=3D 1, y =3D y.'; endif
endfunction

--=20
___________________________________________________________
Sign-up for Ads Free at Mail.com
http://promo.mail.com/adsfreejump.htm