Getting different results when fourier transforming in python vs mathematica?

Question

I am interested in troubleshooting why I get a very different answer when doing fourier transforms of my data in Mathematica vs Python. The version in Mathematica is well behaved, yet all versions of fourier transforms in Python (scipy/numpy) do not seem to work.

Is there any obvious reason why this is the case?

I have attached a somewhat long explanation of my setup, if anyone is interested in the specific situation that I have, which is below:

I have an experiment that works with an electronic signal working at RF, where I send in 10 nanosecond pulses through an AWG, which is amplified with an electrical amplifier. These pulses then drive an optical signal, which has some slower response as a result -- and consequently has a more distorted output signal.

My goal is to be able to send in a corrected input signal to give me the desired, ideal output shape (a gaussian). I'm trying to do this by constructing a response function, by taking a much faster pulse (100 picoseconds), and recording the response of it. From that I Fourier transform this output to get a map in frequency space for the response of my system.

This response function looks like this:

Then, using this response function, I fourier transform my desired signal, divide it by the response function, and then inverse fourier transform my signal back.

When I do this in Mathematica, I get what I want:

Where the red shows a slightly changed signal to compensate for the response function.

Yet when I do the same thing in python, I get the following output:

Which appears to have a phase-scrambled output.

Here is a part of the Python code used:

import numpy as np
from scipy.fft import*

class Fourier():
    def __init__(self, signal, time):
        self.signal = signal
        self.res =time[1]-time[0]
        self.sampling_rate = 1/self.res
        self.period = self.signal.size/self.sampling_rate       
        self.freq_axis = fftshift(fftfreq(self.signal.size, self.res))
        self.fourier = fftshift(fft.self.signal)

    def fouriertrans(self):
        return [self.freq_axis, self.fourier]

    def amplitude(self):
        return abs(self.fourier)/max(abs(self.fourier))

    def phase(self, degree = False):
        return np.angle(self.fourier, deg = degree)

For the inverse Fourier transform, it would be basically the same but then use scipy.fft.ifft, and recalibrate the time axis... And plot the real part of the resulting signal:

np.real(invfourier_signal)

Answer 1

Without seeing your code we can't tell what's wrong. The Python FFT implementation gets used 1000s of times per day so it's safe to assume that this is not the problem.

I recommend standard debugging procedures: Start with simple signals with known answers and work your way up to the point where things break.

You can start with a simple unit impulse, two unit impulses, a sine wave, Gaussian, 1rst order lowpass filter, etc.