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