Update: Another version from
[SciPy-user] Chirp Z transform
by Paul Kienzle and Nadav Horesh has been merged into SciPy, so you should probably just use that. See scipy.signal.CZT, scipy.signal.czt, etc.
This czt copied from:
Re: [Numpy-discussion] zoom FFT with numpy?
Nadav Horesh
Thu, 15 Mar 2007 01:56:03 -0800
and chirpz copied from:
Re: [Numpy-discussion] zoom FFT with numpy?
Stefan van der Walt
Thu, 15 Mar 2007 23:13:24 -0800
| """Chirp z-Transform. | |
| As described in | |
| Rabiner, L.R., R.W. Schafer and C.M. Rader. | |
| The Chirp z-Transform Algorithm. | |
| IEEE Transactions on Audio and Electroacoustics, AU-17(2):86--92, 1969 | |
| """ | |
| import numpy as np | |
| def chirpz(x,A,W,M): | |
| """Compute the chirp z-transform. | |
| The discrete z-transform, | |
| X(z) = \sum_{n=0}^{N-1} x_n z^{-n} | |
| is calculated at M points, | |
| z_k = AW^-k, k = 0,1,...,M-1 | |
| for A and W complex, which gives | |
| X(z_k) = \sum_{n=0}^{N-1} x_n z_k^{-n} | |
| """ | |
| A = np.complex(A) | |
| W = np.complex(W) | |
| if np.issubdtype(np.complex,x.dtype) or np.issubdtype(np.float,x.dtype): | |
| dtype = x.dtype | |
| else: | |
| dtype = float | |
| x = np.asarray(x,dtype=np.complex) | |
| N = x.size | |
| L = int(2**np.ceil(np.log2(M+N-1))) | |
| n = np.arange(N,dtype=float) | |
| y = np.power(A,-n) * np.power(W,n**2 / 2.) * x | |
| Y = np.fft.fft(y,L) | |
| v = np.zeros(L,dtype=np.complex) | |
| v[:M] = np.power(W,-n[:M]**2/2.) | |
| v[L-N+1:] = np.power(W,-n[N-1:0:-1]**2/2.) | |
| V = np.fft.fft(v) | |
| g = np.fft.ifft(V*Y)[:M] | |
| k = np.arange(M) | |
| g *= np.power(W,k**2 / 2.) | |
| return g |
| ## Copyright (C) 2000 Paul Kienzle | |
| ## | |
| ## This program is free software; you can redistribute it and/or modify | |
| ## it under the terms of the GNU General Public License as published by | |
| ## the Free Software Foundation; either version 2 of the License, or | |
| ## (at your option) any later version. | |
| ## | |
| ## This program is distributed in the hope that it will be useful, | |
| ## but WITHOUT ANY WARRANTY; without even the implied warranty of | |
| ## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the | |
| ## GNU General Public License for more details. | |
| ## | |
| ## You should have received a copy of the GNU General Public License | |
| ## along with this program; if not, write to the Free Software | |
| ## Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA | |
| ## usage y=czt(x, m, w, a) | |
| ## | |
| ## Chirp z-transform. Compute the frequency response starting at a and | |
| ## stepping by w for m steps. a is a point in the complex plane, and | |
| ## w is the ratio between points in each step (i.e., radius increases | |
| ## exponentially, and angle increases linearly). | |
| ## | |
| ## To evaluate the frequency response for the range f1 to f2 in a signal | |
| ## with sampling frequency Fs, use the following: | |
| ## m = 32; ## number of points desired | |
| ## w = exp(-2i*pi*(f2-f1)/(m*Fs)); ## freq. step of f2-f1/m | |
| ## a = exp(2i*pi*f1/Fs); ## starting at frequency f1 | |
| ## y = czt(x, m, w, a); | |
| ## | |
| ## If you don't specify them, then the parameters default to a fourier | |
| ## transform: | |
| ## m=length(x), w=exp(2i*pi/m), a=1 | |
| ## Because it is computed with three FFTs, this will be faster than | |
| ## computing the fourier transform directly for large m (which is | |
| ## otherwise the best you can do with fft(x,n) for n prime). | |
| ## TODO: More testing---particularly when m+N-1 approaches a power of 2 | |
| ## TODO: Consider treating w,a as f1,f2 expressed in radians if w is real | |
| function y = czt(x, m, w, a) | |
| if nargin < 1 || nargin > 4, usage("y=czt(x, m, w, a)"); endif | |
| if nargin < 2 || isempty(m), m = length(x); endif | |
| if nargin < 3 || isempty(w), w = exp(2i*pi/m); endif | |
| if nargin < 4 || isempty(a), a = 1; endif | |
| N = length(x); | |
| if (columns(x) == 1) | |
| k = [0:m-1]'; | |
| Nk = [-(N-1):m-2]'; | |
| else | |
| k = [0:m-1]; | |
| Nk = [-(N-1):m-2]; | |
| endif | |
| nfft = 2^nextpow2(min(m,N)+length(Nk)-1); | |
| Wk2 = w.^(-(Nk.^2)/2); | |
| AWk2 = (a.^-k) .* (w.^((k.^2)/2)); | |
| y = ifft(fft(postpad(Wk2,nfft)).*fft(postpad(x,nfft).*postpad(AWk2,nfft))); | |
| y = w.^((k.^2)/2).*y(1+N:m+N); | |
| endfunction |
| import numarray as N | |
| import numarray.fft as F | |
| def czt(x, m=None, w=None, a=1.0): | |
| """ | |
| Copyright (C) 2000 Paul Kienzle | |
| This program is free software; you can redistribute it and/or modify | |
| it under the terms of the GNU General Public License as published by | |
| the Free Software Foundation; either version 2 of the License, or | |
| (at your option) any later version. | |
| This program is distributed in the hope that it will be useful, | |
| but WITHOUT ANY WARRANTY; without even the implied warranty of | |
| MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the | |
| GNU General Public License for more details. | |
| You should have received a copy of the GNU General Public License | |
| along with this program; if not, write to the Free Software | |
| Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 US | |
| usage y=czt(x, m, w, a) | |
| Chirp z-transform. Compute the frequency response starting at a and | |
| stepping by w for m steps. a is a point in the complex plane, and | |
| w is the ratio between points in each step (i.e., radius increases | |
| exponentially, and angle increases linearly). | |
| To evaluate the frequency response for the range f1 to f2 in a signal | |
| with sampling frequency Fs, use the following: | |
| m = 32; ## number of points desired | |
| w = exp(-2i*pi*(f2-f1)/(m*Fs)); ## freq. step of f2-f1/m | |
| a = exp(2i*pi*f1/Fs); ## starting at frequency f1 | |
| y = czt(x, m, w, a); | |
| If you don't specify them, then the parameters default to a Fourier | |
| transform: | |
| m=length(x), w=exp(2i*pi/m), a=1 | |
| Because it is computed with three FFTs, this will be faster than | |
| computing the Fourier transform directly for large m (which is | |
| otherwise the best you can do with fft(x,n) for n prime). | |
| TODO: More testing---particularly when m+N-1 approaches a power of 2 | |
| TODO: Consider treating w,a as f1,f2 expressed in radians if w is real | |
| """ | |
| # Convenience declarations | |
| ifft = F.inverse_fft | |
| fft = F.fft | |
| if m is None: | |
| m = len(x) | |
| if w is None: | |
| w = N.exp(2j*N.pi/m) | |
| n = len(x) | |
| k = N.arange(m, type=N.Float64) | |
| Nk = N.arange(-(n-1), m-1, type=N.Float64) | |
| nfft = next2pow(min(m,n) + len(Nk) -1) | |
| Wk2 = w**(-(Nk**2)/2) | |
| AWk2 = a**(-k) * w**((k**2)/2) | |
| y = ifft(fft(Wk2,nfft) * fft(x*AWk2, nfft)); | |
| y = w**((k**2)/2) * y[n:m+n] | |
| return y | |
| def next2pow(x): | |
| return 2**int(N.ceil(N.log(float(x))/N.log(2.0))) |
| import numarray as N | |
| import numarray.fft as F | |
| def czt(x, m=None, w=None, a=1.0, axis = -1): | |
| """ | |
| Copyright (C) 2000 Paul Kienzle | |
| This program is free software; you can redistribute it and/or modify | |
| it under the terms of the GNU General Public License as published by | |
| the Free Software Foundation; either version 2 of the License, or | |
| (at your option) any later version. | |
| This program is distributed in the hope that it will be useful, | |
| but WITHOUT ANY WARRANTY; without even the implied warranty of | |
| MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the | |
| GNU General Public License for more details. | |
| You should have received a copy of the GNU General Public License | |
| along with this program; if not, write to the Free Software | |
| Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 US | |
| usage y=czt(x, m, w, a) | |
| Chirp z-transform. Compute the frequency response starting at a and | |
| stepping by w for m steps. a is a point in the complex plane, and | |
| w is the ratio between points in each step (i.e., radius increases | |
| exponentially, and angle increases linearly). | |
| To evaluate the frequency response for the range f1 to f2 in a signal | |
| with sampling frequency Fs, use the following: | |
| m = 32; ## number of points desired | |
| w = exp(-2i*pi*(f2-f1)/(m*Fs)); ## freq. step of f2-f1/m | |
| a = exp(2i*pi*f1/Fs); ## starting at frequency f1 | |
| y = czt(x, m, w, a); | |
| If you don't specify them, then the parameters default to a Fourier | |
| transform: | |
| m=length(x), w=exp(2i*pi/m), a=1 | |
| Because it is computed with three FFTs, this will be faster than | |
| computing the Fourier transform directly for large m (which is | |
| otherwise the best you can do with fft(x,n) for n prime). | |
| TODO: More testing---particularly when m+N-1 approaches a power of 2 | |
| TODO: Consider treating w,a as f1,f2 expressed in radians if w is real | |
| """ | |
| # Convenience declarations | |
| ifft = F.inverse_fft | |
| fft = F.fft | |
| do_transpose = (axis != -1) and (x.rank > 1) # transpose data to make it equivalent to axis=-1 | |
| if axis < 0: | |
| axis += x.rank | |
| if do_transpose: | |
| axes = N.arange(x.rank) | |
| axes[[axis, x.rank-1]] = axes[[x.rank-1, axis]] | |
| x = N.transpose(x, axes) | |
| if m is None: | |
| m = x.shape[-1] | |
| if w is None: | |
| w = N.exp(2j*N.pi/m) | |
| n = x.shape[-1] | |
| k = N.arange(m, type=N.Float64) | |
| Nk = N.arange(-(n-1), m-1, type=N.Float64) | |
| nfft = next2pow(min(m,n) + len(Nk) -1) | |
| Wk2 = w**(-(Nk**2)/2) # length = m + len(x) | |
| AWk2 = a**(-k) * w**((k**2)/2) # length = m | |
| y = ifft(fft(Wk2,nfft) * fft(x * N.resize(AWk2, x.shape), nfft)); | |
| y = N.take(y, range(n,m+n), axis=-1) # [n:m+n] | |
| y = N.resize(w**((k**2)/2), y.shape) * y | |
| if do_transpose: | |
| y.transpose(axes) | |
| return y | |
| def next2pow(x): | |
| return 2**int(N.ceil(N.log(float(x))/N.log(2.0))) |