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