FFT of complex signal in Python resulting in "flipped" spectra compared to FFT in Mathematica

user_15953869 0
Viewed 314

I have a complex signal of which I want to take the FFT of. Using Mathematica I got the following result:

(* Some acquisition params *)
fS = 100. 10^6;
time = 10. 10^-6;
NbrSamp = Round[fS * time];

(* Generate signal *)
w0 = 2 \[Pi] 80 10^6;
ti = Subdivide[0., time, NbrSamp];
sig = Cos[w0 ti] + \[ImaginaryI] Sin[w0 ti];

(* Plot result of FFT *)
ListPlot[Abs@Fourier[sig], Joined -> True, PlotRange -> All, DataRange -> {0, fS}]

figure of obtained spectra

However, when performing the same calculations in Python I get a "flipped" spectra compared to Mathematica. In order to get the same spectra I have to flip the output from the FFT. My code reads:

import numpy as np
import matplotlib.pylab as plt
from scipy.fft import fft

# Some acquisition params
fs = 100e6
time = 10e-6
NbrSamp = round(fs*time)

# Generate signal
w0 = 2*np.pi*80e6
ti = np.arange(NbrSamp+1) / fs
sig = np.cos(w0*ti) + 1j*np.sin(w0*ti)

# calc FFT
sigFFT = fft(sig, norm = 'ortho')
freq = np.arange(0, NbrSamp+1) * fs / NbrSamp

# Plot result of FFT
plt.plot(freq, abs(sigFFT), label = "Actual result")
plt.plot(freq, abs(sigFFT[::-1]), label = "Flipped result")
plt.xlim(min(freq), max(freq))
plt.legend()
plt.show()

Results from Python FFT

My question is then why this is and how can I obtain the same results as in Mathematica using Python, without having to flip the data obtained from the FFT? I will have a large number of such signals and therefore I would like to avoid this flip operation.

wolfram-mathematica fft signal-processing python
2 Answers

From the documentation to Mathematica's Fourier function:

  • Other definitions are used in some scientific and technical fields.
  • Different choices of definitions can be specified using the option FourierParameters.
  • With the setting FourierParameters->{a,b}, the discrete Fourier transform computed by Fourier is <some equation with exp(2πib) in it>.
  • Some common choices for {a,b} are {0,1} (default), {-1,1} (data analysis), {1,-1} (signal processing).
  • The setting b = -1 effectively corresponds to conjugating both input and output lists.

So, by default, Mathematica sets b = 0, but setting b = -1 makes its output match Python's, which uses the more common Signal Processing definition.

To get Python to match Mathematica's result, do

sigFFT = np.conj(fft(np.conj(sig), norm = 'ortho'))

Note that this operation is not the same as flipping the result, as the zero frequency component is not affected (the first element of the output array), only the remaining elements are flipped (and conjugated also).

Also note that, if you need the magnitude of the FFT only, it is not necessary to conjugate the output, as abs(conj(x)) == abs(x).

user_73287820

From docs, Mathematica's Fourier Fourier uses a positively-signed (clockwise-spinned) exponent, whereas most Python libs use negatively-signed (as in Wiki); different uses for different applications.

This has the effect of conjugating the input and output: fft_mathematica(x) == conj(fft_python(conj(x), norm='ortho')). To conjugate the output directly (i.e. FFT bins), swap the positives and negatives:

out = fft(x)
out[1:] = out[1:][::-1]
user_101337970
Related