Skip to content

Instantly share code, notes, and snippets.

@endolith
Last active January 7, 2024 11:29
Show Gist options
  • Save endolith/2783807 to your computer and use it in GitHub Desktop.
Save endolith/2783807 to your computer and use it in GitHub Desktop.
Chirp Z-transforms in Python (by Paul Kienzle, Nadav Horesh, Stefan van der Walt)
"""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)))
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment