2

I have an array of complex numbers (waves), which is ordered in the same way as is returned by fft. If I call ifft on this array, it returns an approximation of the original samples.

I would like to implement ifft myself in python. I've found the formula of IFFT. I implemented it, but it looks a bit different from the ifft result. I tried to fix it by looking at the ifft source, but that is a heavily optimized version, and I was not able to find out how that works

This is what I have so far:

import numpy as np
import matplotlib.pyplot as plt

from scipy.fftpack import fft

# create samples
data = np.zeros(1024)
for k in range(0, 1024):
    data[k] = np.sin(k/20 * np.pi*2)

# plot original samples
#plt.plot(data)
#plt.show()

# apply FFT on samples
waves = fft(data)

# plot FFT result
#plt.plot(np.imag(waves),'r')
#plt.plot(np.real(waves),'b')
#plt.ylabel('value')
#plt.xlabel('period')
#plt.show()

# 
res = np.zeros(1024)
for k in range(0, 1024):
    val = 0.0
    for n in range(0, len(waves)-1):
        # https://dsp.stackexchange.com/a/510/25943
        val += waves[n]*np.exp(-1.j * 2*np.pi * n * k / len(waves)) / len(waves)
    res[k] = val.real


#np implementation
res2 = np.fft.ifft(waves)

plt.plot(data, 'b') # original
plt.plot(res,'g') # my implementation
plt.plot(res2,'y') # np implementation
plt.show()

Maybe the zero frequency term and the negative frequency terms has to be handled differently. I am not sure about that, because it is not mentioned in any description of fourier transformation

Iter Ator
  • 8,226
  • 20
  • 73
  • 164

1 Answers1

1

Just two mistakes here: 

for n in range(0, len(waves)-1):

should be 

for n in range(0, len(waves)):

because range does not include its upper limit (this, along with 0-based indexing, makes implementing FFT-type algorithms a bit easier compared to Matlab).

Also,

val += waves[n]*np.exp(-1.j * 2*np.pi * n * k / len(waves))

should be 

val += waves[n]*np.exp(1.j * 2*np.pi * n * k / len(waves))

Sign convention vary; in NumPy, direct transform has -1j and inverse has 1j.

Of course, the whole thing is very inefficient, but you probably wanted to write it out in detail for yourself. NumPy's vectorized operations would be used normally, beginning with 

data = np.zeros(1024)
for k in range(0, 1024):
    data[k] = np.sin(k/20 * np.pi*2)

being replaced by 

data = np.sin(np.arange(1024)/20 * np.pi*2)

and similarly for other loops.

  • How is it possible to do extrapolation correctly? If I just extend the data array and fill it `for k in range(0,2048)` then there is a sudden drop around 1024: https://i.imgur.com/qk5wNnM.png – Iter Ator Oct 29 '17 at 16:46
  • This is a separate, and mathematical question; I suggest asking at [dsp.se], with a clarification of what kind of extrapolation is desired. –  Oct 29 '17 at 17:30