Blog Archive

Friday, October 9, 2015

Python DFT IDFT Studying Note

https://wiki.python.org/moin/NumericAndScientificRecipes

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