2

I have a csv file that I read in (using python 3 on a Jupyter notebook but get the same results from terminal). I am computing the fft via numpy.fft.fft module and am getting the strange result that the fft of the data returns back the original data - i.e a complex vector with real part exactly equal to the (real) input data and the imaginary part identically equal to 0. The code is shown below:

with open('/Users/amacrae/Documents/PMDi/MCT/Jan10/msin287.csv', 'r') as f:
    c = csv.reader(f)
    y = np.array(list(c),dtype=float)
YF = np.fft.fft(y)
print(np.sum(YF.real-y))
print(np.sum(YF.imag))
> 0.0
> 0.0

To ensure that it's not just the data, I plotted the identical data in matlab with the correct results (the data is designed so that the magnitude of the fft is constant over a window in frequency space and has a real iffy.) The corresponding matlab code is:

y = csvread('/Users/amacrae/Documents/PMDi/MCT/Jan10/msin287.csv');
plot(abs(fft(y)))

As far as I can tell the results should be the same in either language ... the real parts of the imported data match in both cases (same length and values) but the fft's do not. The data are quite long - 100,000 samples, but if I create a random 100,000 sample array in python I get a real + imaginary fft. Does anyone have an idea what might be causing this?

Image 1 is in Python. The fft'd data is precisely the input data

Image 2 is the same data, but fft'd in matlab

Will Vousden
  • 32,488
  • 9
  • 84
  • 95
Orko
  • 142
  • 1
  • 7

2 Answers2

2

Reagarding your code

print(np.sum(YF.real-y))
print(np.sum(YF.imag))

The first part is true because of Parseval's theorem The second is true since the first half of the spectrum and the second half of the spectrum are conjugate.

Try comparing the absolute value in Python with your Matlab version. Absolute value differs from real. In both cases consider defining the length of the transform i.e. 1024, 2048 (since e^(j*1e-5) may cause problem)

Will Vousden
  • 32,488
  • 9
  • 84
  • 95
Radu Ionescu
  • 3,462
  • 5
  • 24
  • 43
  • Hi, thx for the reply. Actually this is not true. First, the first and second parts of y are not complex conjugates so fft(y) needn't be real. In fact in this case, the spectrum is complex by construction with phases necessary to produce a real sequence. Second, Parseval's theorem does not guarantee this. What it does say is that the powers of y and YF are equal: sum(abs(y)**2) = sum(abs(Y)**2). For example, with y = [1,2,3], YF = 1.5*[4, -1+1i/sqrt(3), -1-1i/sqrt(3)]. Here, sum(YF.real-y) = -3, not 0. So something else must be going on here.... – Orko Jan 12 '16 at 16:40
1

Thanks for the input Radu. In the end it turned out to be something else (see comment above.)

The issue was that the method of loading the file produced an array of arrays, i.e. a column vector. When I called fft, it just returned the original column vector back as it performed the fft on each row individually.

with open('/Users/amacrae/Documents/PMDi/MCT/Jan10/msin287.csv', 'r') as f:
    c = csv.reader(f)
    y = np.array(list(c),dtype=float)
    # y = [[y0],[y1],[y2],...]] fft(y) = [[y0+0.j],[y1+0.j],[y2+0.j],...]]

The solution was to either flatten the array:

with open('/Users/amacrae/Documents/PMDi/MCT/Jan10/msin287.csv', 'r') as f:
    c = csv.reader(f)
    y = np.array(list(c),dtype=float)
    y2 = flatten(y)
    # produces [y0,y1,y2,....]

or use a different method for reading the file:

y2 = np.loadtxt('/Users/amacrae/Documents/PMDi/MCT/Jan10/msin287.csv')
# also produces [y0,y1,y2,....]

Both produced the correct FFT.

Orko
  • 142
  • 1
  • 7