Note:
cmath — Mathematical functions for complex numbers
"""
def fft(x, sign=-1):
from cmath import pi, exp
N = len(x)
W = [exp(sign * 2j * pi * i / N)
for i in range(N)] # exp(-j...) is default
x = bitrev(x)
m = 2
while m <= N:
for s in range(0, N, m):
for i in range(m/2):
n = i * N / m
a, b = s + i, s + i + m/2
x[a], x[b] = x[a] + W[n % N] * x[b], x[a] - W[n % N] * x[b]
m *= 2
return x
"""
Inverse FFT with normalization by N, so that x == ifft(fft(x)) within
round-off errors.
"""
def ifft(X):
N, x = len(X), fft(X, sign=1) # e^{j2\pi/N}
for i in range(N):
x[i] /= float(N)
return x
"""
DFT using discrete summation
X(n) = \sum_k W^{nk} x(k), W = e^{-j2\pi/N}
where N need not be power of 2. The choice of e^{-j2\pi/N} or
e^{j2\pi/N} is made by "sign=-1" or "sign=1" respectively.
"""
def dft(x, sign=-1):
from cmath import pi, exp
N = len(x)
W = [exp(sign * 2j * pi * i / N)
for i in range(N)] # exp(-j...) is default
X = [sum(W[n * k % N] * x[k] for k in range(N))
for n in range(N)]
return X
"""
Inverse DFT with normalization by N, so that x == idft(dft(x)) within
round-off errors.
"""
def idft(X):
N, x = len(X), dft(X, sign=1) # e^{j2\pi/N}
for i in range(N):
x[i] /= float(N)
return x
"""
No comments:
Post a Comment