Making a function to take the inverse fourier transform

Question

I am trying to take the inverse fourier transform by making my own function. This is the function to take the Fourier transform of my time series which appears to work fine.

def DFT(x, frequencies):
    N1 = x.size
    k = frequencies.size
    t = np.arange(N1)
    fourier = np.zeros(k)
    for i in range(0,k):
        fourier[i] = np.dot(x,(1/k)*np.exp(-2j*np.pi*frequencies[i]*t))
    return fourier

This is my original signal (just a sine wave):

N2 = 1*10**6
time = np.arange(0,N2,1000)
lam = .1*10**6
f = 1/lam
freq = np.arange(0,.05,.0001)


signal = np.sin(2*np.pi*f*time)

And the power spectrum is plotted using my DFT (fourier function):

plt.plot(freq, np.abs(DFT(signal,freq))**2)
plt.xlabel('Frequency')
plt.title('Spectrum of Sine Wave')
plt.grid()
plt.show()

but when I try to apply my function for the inverse fourier transform, I am not getting my original sine wave back:

def IFT(fft, frequencies):
    N = fft.size
    k = frequencies.size
    n = np.arange(N)
    inverse_fourier = np.zeros(k)
    for i in range(0,k):
        inverse_fourier[i] = np.dot(fft,np.exp((-2j*np.pi*frequencies[i]*n)/N)) #[None,:]
    return inverse_fourier

What is wrong with my function? I get no errors, but the returned signal is totally wrong.

Answer 1

Running you code you should get the following warning:

ComplexWarning: Casting complex values to real discards the imaginary part
  fourier[i] = np.dot(x,(1/k)*np.exp(-2j*np.pi*frequencies[i]*t))

Since the resulting Fourier transform should be complex valued, this warning should be reasons for concerns. To get rid of this warning you may initialize fourier like so:

fourier = np.zeros(k, dtype=complex)

Also the formula for Discrete Fourier Transform includes summations over frequencies covering the complete [0,1) range. To get a 1000-point DFT (as you had in your code) you'd then have to use

freq = np.arange(0,1,.001)

This will result in a spectrum that includes 2 spikes: one at the expected frequency, and another symmetric one above the Nyquist frequency. It is common to discard the results above the Nyquist frequency when plotting the spectrum of real-valued signals (but use the full spectrum into your IFT function).

Finally, as GrimTrigger pointed out:

your inverse the exponent should be positive (2j instead of -2j) and drop the /N

Answer 2

In your inverse the exponent should be positive (2j instead of -2j) and drop the /N, which gives (added plots for demonstration):

import numpy as np
import matplotlib.pyplot as plt

def DFT(x, frequencies):
    N1 = x.size
    k = frequencies.size
    t = np.arange(N1)
    fourier = np.zeros(k)
    for i in range(0,k):
        fourier[i] = np.dot(x, (1/k)*np.exp(-2j*np.pi*frequencies[i]*t))
    return fourier

def IFT(fft, frequencies):
    N = fft.size
    k = frequencies.size
    n = np.arange(N)
    inverse_fourier = np.zeros(k)
    for i in range(0,k):
        inverse_fourier[i] = np.dot(fft, np.exp((2j*np.pi*frequencies[i]*n))) #[None,:]
    return inverse_fourier

N2 = 1*10**6
time = np.arange(0,N2,2000)
lam = .1*10**6
f = 1/lam
freq = np.arange(0,.05,.0001)

signal = np.sin(2*np.pi*f*time)

plt.plot(time, signal)
plt.xlabel('Time')
plt.title('Sine Wave')
plt.grid()
plt.show()

dft = DFT(signal, freq)

plt.plot(freq, np.abs(dft)**2)
plt.xlabel('Frequency')
plt.title('Spectrum of Sine Wave')
plt.grid()
plt.show()

plt.plot(time, IFT(dft, freq))
plt.xlabel('Time')
plt.title('Sine Wave')
plt.grid()
plt.show()

which gives (first sin graph omitted):

and