Mercurial > forge
comparison main/signal/czt.m @ 0:6b33357c7561 octave-forge
Initial revision
| author | pkienzle |
|---|---|
| date | Wed, 10 Oct 2001 19:54:49 +0000 |
| parents | |
| children | 6cd6668c225b |
comparison
equal
deleted
inserted
replaced
| -1:000000000000 | 0:6b33357c7561 |
|---|---|
| 1 ## Copyright (C) 2000 Paul Kienzle | |
| 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 2 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; if not, write to the Free Software | |
| 15 ## Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA | |
| 16 | |
| 17 ## usage y=czt(x, m, w, a) | |
| 18 ## | |
| 19 ## Chirp z-transform. Compute the frequency response starting at a and | |
| 20 ## stepping by w for m steps. a is a point in the complex plane, and | |
| 21 ## w is the ratio between points in each step (i.e., radius increases | |
| 22 ## exponentially, and angle increases linearly). | |
| 23 ## | |
| 24 ## To evaluate the frequency response for the range f1 to f2 in a signal | |
| 25 ## with sampling frequency Fs, use the following: | |
| 26 ## m = 32; ## number of points desired | |
| 27 ## w = exp(-2i*pi*(f2-f1)/(m*Fs)); ## freq. step of f2-f1/m | |
| 28 ## a = exp(2i*pi*f1/Fs); ## starting at frequency f1 | |
| 29 ## y = czt(x, m, w, a); | |
| 30 ## | |
| 31 ## If you don't specify them, then the parameters default to a fourier | |
| 32 ## transform: | |
| 33 ## m=length(x), w=exp(2i*pi/m), a=1 | |
| 34 ## Because it is computed with three FFTs, this will be faster than | |
| 35 ## computing the fourier transform directly for large m (which is | |
| 36 ## otherwise the best you can do with fft(x,n) for n prime). | |
| 37 | |
| 38 ## TODO: More testing---particularly when m+N-1 approaches a power of 2 | |
| 39 ## TODO: Consider treating w,a as f1,f2 expressed in radians if w is real | |
| 40 function y = czt(x, m, w, a) | |
| 41 if nargin < 1 || nargin > 4, usage("y=czt(x, m, w, a)"); endif | |
| 42 if nargin < 2 || isempty(m), m = length(x); endif | |
| 43 if nargin < 3 || isempty(w), w = exp(2i*pi/m); endif | |
| 44 if nargin < 4 || isempty(a), a = 1; endif | |
| 45 | |
| 46 N = length(x); | |
| 47 if (columns(x) == 1) | |
| 48 k = [0:m-1]'; | |
| 49 Nk = [-(N-1):m-2]'; | |
| 50 else | |
| 51 k = [0:m-1]; | |
| 52 Nk = [-(N-1):m-2]; | |
| 53 endif | |
| 54 nfft = 2^nextpow2(min(m,N)+length(Nk)-1); | |
| 55 Wk2 = w.^(-(Nk.^2)/2); | |
| 56 AWk2 = (a.^-k) .* (w.^((k.^2)/2)); | |
| 57 y = ifft(fft(postpad(Wk2,nfft)).*fft(postpad(x,nfft).*postpad(AWk2,nfft))); | |
| 58 y = w.^((k.^2)/2).*y(1+N:m+N); | |
| 59 endfunction |