Reconstruct original signal with FFT in python

Question

I have an analitically generated spectrum, where x axis represents angular frequency, y represents intensity. The spectrum is centered around some frequency value, which is often called central frequency of the signal (blue graph on the picture). I want to perform IFFT on the data to time domain, cut its useful part with a gaussian curve, then FFT back to the original domain . My problem is that after IFFT(FFT(signal)) the central frequency is lost: I get back the spectrum by shape, but it's always centered around 0 (orange graph). Currently my solution to this is quite bad: I cache the original x axis and I restore it upon FFT calls. This obviously has many downsides, and I want to improve it. Below I included a small demo which demonstrates the problem. My question is: can this be solved in a more elegant way? Is there a way central frequency is not lost during the process?

import numpy as np
from scipy.fftpack import fft, ifft, fftshift, fftfreq
from scipy.interpolate import interp1d
import matplotlib.pyplot as plt

C_LIGHT = 299.793

def generate_data(start, stop, center, delay, GD=0, resolution=0.1):
    window = 8 * np.log(2) / 50
    lamend = 2 * np.pi * C_LIGHT / start
    lamstart = 2 * np.pi * C_LIGHT/stop
    lam = np.arange(lamstart, lamend + resolution, resolution) 
    omega = 2 * np.pi * C_LIGHT / lam 
    relom = omega - center
    _i = np.exp(-(relom) ** 2 / window)
    i = 2 * _i + 2 * np.cos(relom * GD + (omega * delay)) * np.sqrt(_i * _i)
    return omega, i


if __name__ == '__main__':

    # Generate data
    x, y = generate_data(1, 3, 2, 800, GD=0)

    # Linearly interpolate to be evenly spaced
    xs = np.linspace(x[0], x[-1], len(x))
    intp = interp1d(x, y, kind='linear')
    ys = intp(xs)
    x, y = xs, ys
    plt.plot(x, y, label='original')

    # IFFT 
    xt = fftfreq(len(x), d=(x[0]-x[1])/(2*np.pi))
    yt = ifft(y)
    # plt.plot(xt, np.abs(yt))

    # FFT back
    xf = fftshift(fftfreq(len(xt), d=(xt[0]-xt[1])/(2*np.pi)))
    yf = fft(yt)
    plt.plot(xf, np.abs(yf), label='after transforms')
    plt.legend()
    plt.grid()
    plt.show()

Answer 1

I think that fftfreq does not do what you think it does. The xf for fft(ifft(y) is identical to x, you should not try to re-compute it. The x-axis doesn't change when going to another domain and then back again.

Also, do note that fftfreq returns the coordinates in the frequency domain for the discrete Fourier transform of a signal of a given length and with a given sample spacing. It does not do the reverse, you cannot use it to determine the coordinates in the spatial domain after applying the inverse discrete Fourier transform. (The spacing it returns is valid, but the set of coordinates is not.)

    plt.plot(x, y, label='original')

    # IFFT 
    yt = ifft(y)
    # plt.plot(np.abs(yt))

    # FFT back
    yf = fft(yt)
    plt.plot(x, np.real(yf), label='after transforms')
    plt.legend()
    plt.grid()
    plt.show()

---

Another problem with your code is that ifft(y) assumes a fixed set of values along the x-axis. Your x does not match this. Thus, the spatial-domain signal you obtain is not meaningful.

Running your code, I see that x runs from 3.0 to 1.0 in steps of 0.0004777. You will have to augment your data so that the values run from 0.0 to 6.0, with the region (3.0, 6.0) being the conjugate symmetric copy of the region (0.0, 3.0). This region corresponds to the negative frequencies, according to the periodicity of the frequency domain (F[n]==F[n+N], with N the number of samples). Fill the region (0.0, 1.0) with zeros.

Given this standardized x-axis in the frequency domain, xf = fftfreq(len(xt), d=(xt[1]-xt[0])) should reconstruct the x-axis. But you need to compute xt appropriately: xt = np.linspace(0, 1/(x[1]-x[0]), len(x), endpoint=False) (with x the standardized DFT frequency axis).