numpy.fft 实现 czt (Chirp Z-transform)

动机

如果对L2(R)L^2(\R)L2(R)上做Fourier变换，直接用离散FFT是不行的。需要用CZT。用于数值计算的numpy没有提供CZT，需要重新实现。本文用FFT实现CZT。

f^(ω)=∫Rf(x)e−2πiωxdx∼∫[a,b]f(x)e−2πiωxdx,−a,b≫1∼b−aN(∑l=0N−1f(lN)e−2πiω(a+(b−a)l/N)+f(b)e−2πiωb−f(a)e−2πiωa2)=b−aNe−2πiωa(∑l=0N−1f(lN)e−2πi(b−a)ωl/N+f(b)e2πiω(a−b)−f(a)2)⟹f^(c+kd−cM)=b−aNe−2πia(c+kd−cM)(∑l=0N−1f(lN)e−2πib−aNcle−2πi(b−a)(d−c)MNkl+...)=b−aNe−2πia(c+kd−cM)(CZT({fl},M,e−2πi(b−a)(d−c)MN,e2πib−aNc)+...)\hat{f}(\omega)=\int_\R f(x)e^{-2\pi i\omega x} \mathrm{d}x\\ \sim \int_{[a, b]} f(x)e^{-2\pi i\omega x} \mathrm{d}x, -a, b\gg 1 \\ \sim \frac{b-a}{N}(\sum_{l=0}^{N-1} f(\frac{l}{N})e^{-2\pi i \omega (a+(b-a)l/N)} + \frac{f(b)e^{-2\pi i\omega b}-f(a)e^{-2\pi i\omega a}}{2})\\ = \frac{b-a}{N}e^{-2\pi i \omega a}(\sum_{l=0}^{N-1} f(\frac{l}{N})e^{-2\pi i (b-a) \omega l/N} + \frac{f(b)e^{2\pi i\omega (a-b)}-f(a)}{2})\\ \Longrightarrow\\ \hat{f}(c+k\frac{d-c}{M})=\frac{b-a}{N}e^{-2\pi i a(c+k\frac{d-c}{M})}(\sum_{l=0}^{N-1} f(\frac{l}{N})e^{-2\pi i \frac{b-a}{N} cl}e^{-2\pi i \frac{(b-a)(d-c)}{MN}k l} + ...)\\ =\frac{b-a}{N}e^{-2\pi i a(c+k\frac{d-c}{M})}(CZT(\{f_l\}, M, e^{-2\pi i \frac{(b-a)(d-c)}{MN}}, e^{2\pi i \frac{b-a}{N} c})+...) f^(ω)=∫Rf(x)e−2πiωxdx∼∫[a,b]f(x)e−2πiωxdx,−a,b≫1∼Nb−a(l=0∑N−1f(Nl)e−2πiω(a+(b−a)l/N)+2f(b)e−2πiωb−f(a)e−2πiωa)=Nb−ae−2πiωa(l=0∑N−1f(Nl)e−2πi(b−a)ωl/N+2f(b)e2πiω(a−b)−f(a))⟹f^(c+kMd−c)=Nb−ae−2πia(c+kMd−c)(l=0∑N−1f(Nl)e−2πiNb−acle−2πiMN(b−a)(d−c)kl+...)=Nb−ae−2πia(c+kMd−c)(CZT({fl},M,e−2πiMN(b−a)(d−c),e2πiNb−ac)+...)

原理

基于Bluestein算法
参考 Chirp Z-transform 及其中引文

源码

def czt(x, m=None, w=None, a=1):"""chirp z transform: (Based on FFT)z = A * W.^(-(0:M-1)), g_i = sum_jz^{-ij}x_jInputs m, w and a must be scalars.Arguments:x {array} -- the original signalKeyword Arguments:m {number} -- (default: {None})w {complex number} -- base of the series (default: {np.exp(-2j * np.pi / m)})a {number} -- (default: {1})Returns:array -- czt of x"""n = x.shape[0]if k is None:m = nif w is None:w = np.exp(-2j * np.pi / m)if not isinstance(x, np.ndarray):x = np.asarray(x, dtype=np.complex)# Length for power-of-two fft.nfft = 2 ** nextpow2(n+k-1)# Premultiply data.k = np.arange(-n+1, max((m, n))) ** 2 / 2ww = w ** k   # Chirp filter is 1./ww# nn = n-1 : 2*n-1aa = a ** (-np.arange(0, n)) * ww[n-1:2*n-1]y = x * aa# Fast convolution via FFT.fy = np.fft.fft(y, nfft) * np.fft.fft(1 / ww[:m+n-1], nfft)g = np.fft.ifft(fy)# Finally multiply Chirp filter.return g[n-1:n+m-1] * ww[n-1:n+m-1]

Fourier 变换实现

有了czt就可以实现FT了

def cft(f, trange=(0,1), Nt=None, wrange=(-10,10), Nw=101, padding='c'):"""continuous Fourier transformAssertion: wsize * tsize <= Nt;Reference: https://en.wikipedia.org/wiki/Chirp_Z-transformArguments:f {function|array} -- function or values of functionsKeyword Arguments:trange {tuple} -- the range of time (default: {(0,1)})Nt {[type]} -- number of samples on time-domain (default: {None})wrange {tuple} -- the range of frequency (default: {(-10,10)})Nw {number} -- number of samples on frequency-domain (default: {101})padding {str} -- Used in feature (default: {'c'})Returns:tuple -- function values and frequencies"""a, b = trangec, d = wrangedt = (b-a)/ Nt; dw = (d-c)/ Nwif callable(f):if Nt is None:Nt = 100f = f(np.linspace(a, b, Nt))else:if Nt is None:Nt = f.shape[0]else:f = f[:Nt]F = czt(f, Nw-1, np.exp(-2j*np.pi*dt*dw), np.exp(2j*np.pi*dt*c))ws = np.linspace(c, d, Nw)if padding=='c':p = f[-1]elif padding=='c1':p = 2*f[-1]-f[-2]F = dt * np.exp(-2j*np.pi*a*ws[:-1]) * (F + (p * np.exp(2j*np.pi*(a-b) * ws[:-1]) - f[0])/2)return np.concatenate([F, [F[-1]]]), ws