0
$\begingroup$

I have an analitically generated spectrum $I(\omega)$, 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. I want to perform IFFT on the data to time domain $I(\tau)$, cut its useful part with a gaussian curve, then FFT back to the original domain $I(\omega)$. My problem is that after $F(F^{-1}(I(\omega)))$ the central frequency is lost. I get back the spectrum by shape, but it's always centered around 0. Spectrum 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? Am I missing something?

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()
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.