|
0
|
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 |
