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