Why is every second coefficient of numpy's fft the inverse of what it should be?

Question

This is how I tried to get the DFT of the unit pulse using numpy (the plot shows the unit pulse):

%matplotlib inline
import matplotlib.pyplot as plt
import numpy as np

def plot_complex(space, arr):
    plt.figure()
    plt.plot(space, arr.real, label="real")
    plt.plot(space, arr.imag, label="imag")
    plt.legend(loc='upper left')

f = lambda x: 1 if abs(x) < 0.5 else 0
lo = -10
hi = 10
num_samples = 1000
sample_freq = num_samples/(hi-lo)
linspace = np.linspace(lo, hi, num_samples)
array = np.vectorize(f)(linspace)

coeff = np.fft.fft(array, num_samples)

# we need to shift the coefficients because otherwise negative frequencies are in the right half of the array
coeff_shifted = np.fft.fftshift(coeff)

plot_complex(linspace, array)
plt.title("The unit pulse")

Unit pulse plot

This appears to be working, because using the inverse fft recovers the original unit pulse:

icoeff = np.fft.ifft(coeff)

plot_complex(linspace, icoeff)
plt.title("The recovered signal")

Recovered pulse plot

However, when looking at the signal, it doesn't look like the sinc function I would expect from looking at the continuous fourier transform:

freqspace = np.vectorize(lambda x: x * sample_freq)(np.fft.fftshift(np.fft.fftfreq(1000)))
plot_complex(freqspace, coeff_shifted)
plt.title("DFT coefficients")

plot_complex(freqspace[450:550], coeff_shifted[450:550])
plt.title("Zoomed in")

Plot of the coefficients

Zoomed in plot

It only looks like the sinc function once every second coefficient is multiplied by -1:

# multiplies every second number in the array by -1
def flip_seconds(coeff):
    return np.array([(1 if i%2 == 0 else -1) * s for (i,s) in enumerate(coeff)])

plot_complex(freqspace, flip_seconds(coeff_shifted))

Plot of the altered coefficients

Why is this?

It sounds like your inverse FFT doesn't have the phase information. I'm not sure on the numpy specifics, but take a look at this answer https://dsp.stackexchange.com/a/32366/7786 — Peter Gibson, Aug 02 '18 at 03:47
@PeterGibson I'd be rather surprised if I lost phase information somehow - my code seems to be able to recover any function (real or complex) essentially perfectly — nnotm, Aug 02 '18 at 03:54
You aren't accounting for the phase properly. Plot the *magnitude* of the complex Fourier coefficients, and you'll see the absolute value of the sinc function. — Warren Weckesser, Aug 02 '18 at 03:55
@WarrenWeckesser okay, I suppose that makes sense, but why does the phase flip after every coefficient rather than being continuous, like in the continuous FT? — nnotm, Aug 02 '18 at 04:03
To produce the case where the imaginary parts of the Fourier coefficients are all exactly zero, you have position the pulse properly in the input sequence. For example, if `x = np.array([1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1])`, and you compute `y = np.fft.fft(x)`, then `y.imag` is all 0 (and `y.real` will be the sinc function). (Use `fftshift` as desired to center the plot where you like.) — Warren Weckesser, Aug 02 '18 at 04:18
(When I say the "sinc function", I really mean the "periodic sinc function", aka the [Dirichlet kernel](https://en.wikipedia.org/wiki/Dirichlet_kernel).) — Warren Weckesser, Aug 02 '18 at 04:27
@WarrenWeckesser Ah, I hadn't realized that the pulse needs to span over the edges rather than be in the middle! It makes sense in retrospect, given that this is how the frequency space works. It works if I call np.fft.fftshift on the array after creating it. I suppose if you want to write this as an answer, I'll accept it. — nnotm, Aug 02 '18 at 04:28
@nnotm I see that you got to the answer while I was I was writing up my reply :). — Bi Rico, Aug 02 '18 at 04:34

score 2 · Accepted Answer · answered Aug 02 '18 at 04:33

This is a bit handwavy, maybe someone else wants to flush it out with the math :) but basically you've windowed your "pulse domain" to be [-X/2, X/2] while the fft expects it to be windowed [0, X]. The difference is a shift in the "pulse domain", which results in a phase shift in the frequency domain. Because you've shifted by exactly half the window, the phase shift is exactly exp(-pi * f * i / sampling_freq) (or something like that) so it shows up as every other item being multiplied as exp(-pi * i). You can fix this by shifting your "pulse space" before applying the fft.

import matplotlib.pyplot as plt
import numpy as np

def plot_complex(space, arr):
    plt.figure()
    plt.plot(space, arr.real, label="real")
    plt.plot(space, arr.imag, label="imag")
    plt.legend(loc='upper left')

lo = -10
hi = 10
num_samples = 1000
sample_freq = num_samples/(hi-lo)
linspace = np.linspace(lo, hi, num_samples)
array = abs(linspace) < .5
array_shifted = np.fft.fftshift(array)

coeff = np.fft.fft(array_shifted, num_samples)

# we need to shift the coefficients because otherwise negative frequencies are in the right half of the array
coeff_shifted = np.fft.fftshift(coeff)

plot_complex(linspace, array_shifted)
plt.title("The unit pulse (shifted)")

icoeff = np.ftt.ifftshift(np.fft.ifft(coeff))

plot_complex(linspace, icoeff)
plt.title("The recovered signal")

freqspace = np.fft.fftshift(np.fft.fftfreq(1000)) * sample_freq
plot_complex(freqspace, coeff_shifted)
plt.title("DFT coefficients")

plot_complex(freqspace[450:550], coeff_shifted[450:550])
plt.title("Zoomed in")


def flip_seconds(coeff):
    return np.array([(1 if i%2 == 0 else -1) * s for (i,s) in enumerate(coeff)])

plot_complex(freqspace, flip_seconds(coeff_shifted))
plt.show()

You need to use `np.fft.ifftshift` to shift the signal before computing the FFT. `fftshift` shifts the origin from the leftmost bin to the middle, `ifftshift` shifts it from the middle to the leftmost bin. These two operations are not the same in case the signal is odd in length. — Cris Luengo, Aug 02 '18 at 05:10

Why is every second coefficient of numpy's fft the inverse of what it should be?

1 Answers1