octave_packages/signal-1.1.3/czt.m

   1 ## Copyright (C) 2004 Daniel Gunyan
   2 ##
   3 ## This program is free software; you can redistribute it and/or modify it under
   4 ## the terms of the GNU General Public License as published by the Free Software
   5 ## Foundation; either version 3 of the License, or (at your option) any later
   6 ## version.
   7 ##
   8 ## This program is distributed in the hope that it will be useful, but WITHOUT
   9 ## ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
  10 ## FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
  11 ## details.
  12 ##
  13 ## You should have received a copy of the GNU General Public License along with
  14 ## this program; if not, see <http://www.gnu.org/licenses/>.
  15
  16 ## usage y=czt(x, m, w, a)
  17 ##
  18 ## Chirp z-transform.  Compute the frequency response starting at a and
  19 ## stepping by w for m steps.  a is a point in the complex plane, and
  20 ## w is the ratio between points in each step (i.e., radius increases
  21 ## exponentially, and angle increases linearly).
  22 ##
  23 ## To evaluate the frequency response for the range f1 to f2 in a signal
  24 ## with sampling frequency Fs, use the following:
  25 ##     m = 32;                               ## number of points desired
  26 ##     w = exp(-j*2*pi*(f2-f1)/((m-1)*Fs));  ## freq. step of f2-f1/m
  27 ##     a = exp(j*2*pi*f1/Fs);                ## starting at frequency f1
  28 ##     y = czt(x, m, w, a);
  29 ##
  30 ## If you don't specify them, then the parameters default to a fourier
  31 ## transform:
  32 ##     m=length(x), w=exp(-j*2*pi/m), a=1
  33 ##
  34 ## If x is a matrix, the transform will be performed column-by-column.
  35
  36 ## Algorithm (based on Oppenheim and Schafer, "Discrete-Time Signal
  37 ## Processing", pp. 623-628):
  38 ##   make chirp of length -N+1 to max(N-1,M-1)
  39 ##     chirp => w^([-N+1:max(N-1,M-1)]^2/2)
  40 ##   multiply x by chirped a and by N-elements of chirp, and call it g
  41 ##   convolve g with inverse chirp, and call it gg
  42 ##     pad ffts so that multiplication works
  43 ##     ifft(fft(g)*fft(1/chirp))
  44 ##   multiply gg by M-elements of chirp and call it done
  45
  46 function y = czt(x, m, w, a)
  47   if nargin < 1 || nargin > 4, print_usage; endif
  48
  49   [row, col] = size(x);
  50   if row == 1, x = x(:); col = 1; endif
  51
  52   if nargin < 2 || isempty(m), m = length(x(:,1)); endif
  53   if length(m) > 1, error("czt: m must be a single element\n"); endif
  54   if nargin < 3 || isempty(w), w = exp(-2*j*pi/m); endif
  55   if nargin < 4 || isempty(a), a = 1; endif
  56   if length(w) > 1, error("czt: w must be a single element\n"); endif
  57   if length(a) > 1, error("czt: a must be a single element\n"); endif
  58
  59   ## indexing to make the statements a little more compact
  60   n = length(x(:,1));
  61   N = [0:n-1]'+n;
  62   NM = [-(n-1):(m-1)]'+n;
  63   M = [0:m-1]'+n;
  64
  65   nfft = 2^nextpow2(n+m-1); # fft pad
  66   W2 = w.^(([-(n-1):max(m-1,n-1)]'.^2)/2); # chirp
  67
  68   for idx = 1:col
  69     fg = fft(x(:,idx).*(a.^-(N-n)).*W2(N), nfft);
  70     fw = fft(1./W2(NM), nfft);
  71     gg = ifft(fg.*fw, nfft);
  72
  73     y(:,idx) = gg(M).*W2(M);
  74   endfor
  75
  76   if row == 1, y = y.'; endif
  77 endfunction
  78
  79 %!shared x
  80 %! x = [1,2,4,1,2,3,5,2,3,5,6,7,8,4,3,6,3,2,5,1];
  81 %!assert(fft(x),czt(x),10000*eps);
  82 %!assert(fft(x'),czt(x'),10000*eps);
  83 %!assert(fft([x',x']),czt([x',x']),10000*eps);