# endolith / chirpz.py Last active March 11, 2015

### SSH clone URL

You can clone with HTTPS or SSH.

Chirp Z-transforms in Python (by Paul Kienzle, Nadav Horesh, Stefan van der Walt)

View chirpz.py
 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 """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 
View chirpz.py
 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 ## 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 
View chirpz.py
 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 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))) 
View chirpz.py
 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 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))) 
View chirpz.py
to join this conversation on GitHub. Already have an account? Sign in to comment
Something went wrong with that request. Please try again.