0
$\begingroup$

I edited the post to be more specific, see below.

I have a dataset, on which I need to perform and IFFT, cut the valueable part of it (by multiplying with a gaussian curve), then FFT back. First it's in angular frequency domain, so an IFFT leads to time domain. Then FFT-ing back should lead to angular frequency again, but I can't seem to find a solution how to get back the original domain.

For the transforms I'm using np.fft.fftfreq the following way:

# x is in ang. frequency domain, that's the reason for the 2*np.pi division
t = np.fft.fftfreq(len(x), d=(x[1]-x[0])/(2*np.pi))

However doing

x = np.fft.fftfreq(len(t), d=2*np.pi*(t[1]-t[0]))

completely not giving me back the original x values. Is that something I'm misunderstanding?

The question can be asked generalized, for example:

import numpy as np

x = np.arange(100)

xx = np.fft.fftfreq(len(x), d = x[1]-x[0])
# how to get back the original x after fftfreq?

I've tried to use a temporal variable where I store the original x values, but it's not too elegant. I'm looking for some kind of inverse of fftfreq, and in general the possible best solution for that problem. Thank you.

EDIT: I provide there a minimal reproducible example. I have a dataset which has angular frequency on x axis and intensity on the y. I want to perfrom IFFT to change to time domain. Unfortunately the x values are not evenly spaced, so a (linear) interpolation is needed first before IFFT. Then in time domain the transform looks like this: See this picture

The next step is to cut one of the symmetrical spikes with a gaussian curve, then FFT back to angular frequency domain (the same where we started). My problem is when I transfrom the x-axis for the IFFT (which I think is correct), I can't get back into the original angular frequency domain. Here is the code, which includes the generator for the dataset too.

import numpy as np
import matplotlib.pyplot as plt 
import scipy
from scipy.interpolate import interp1d

C_LIGHT = 299.792

# for easier case, this is zero, so it can be ignored.
def _disp(x, GD=0, GDD=0, TOD=0, FOD=0, QOD=0):
    return x*GD+(GDD/2)*x**2+(TOD/6)*x**3+(FOD/24)*x**4+(QOD/120)*x**5

# the generator to make sample datasets
def generator(start, stop, center, delay, GD=0, GDD=0, TOD=0, FOD=0, QOD=0, resolution=0.1, pulse_duration=15, chirp=0):
    window = (np.sqrt(1+chirp**2)*8*np.log(2))/(pulse_duration**2)
    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_r = np.exp(-(relom)**2/(window))
    i_s = np.exp(-(relom)**2/(window))

    i = i_r + i_s + 2*np.sqrt(i_r*i_s)*np.cos(_disp(relom, GD=GD, GDD=GDD, TOD=TOD, FOD=FOD, QOD=QOD)+delay*omega)
    #since the _disp polynomial is set to be zero, it's just cos(delay*omega)
    return omega, i


def interpol(x,y):
    ''' Simple linear interpolation '''
    xs = np.linspace(x[0], x[-1], len(x))
    intp = interp1d(x, y, kind='linear', fill_value = 'extrapolate')
    ys = intp(xs)
    return xs, ys


def ifft_method(initSpectrumX, initSpectrumY, interpolate=True):
    if len(initSpectrumY) > 0 and len(initSpectrumX) > 0:
        Ydata = initSpectrumY
        Xdata = initSpectrumX
    else:
        raise ValueError
    N = len(Xdata)
    if interpolate:
        Xdata, Ydata = interpol(Xdata, Ydata)
        # the (2*np.pi) division is because we have angular frequency, not frequency
        xf = np.fft.fftfreq(N, d=(Xdata[1]-Xdata[0])/(2*np.pi)) * N * Xdata[-1]/(N-1)
        yf = np.fft.ifft(Ydata)
    else:
        pass # some irrelevant code there
    return xf, yf

def fft_method(initSpectrumX ,initSpectrumY):
    if len(initSpectrumY) > 0 and len(initSpectrumX) > 0:
        Ydata = initSpectrumY
        Xdata = initSpectrumX
    else:
        raise ValueError
    yf = np.fft.fft(Ydata)
    xf = np.fft.fftfreq(len(Xdata), d=(Xdata[1]-Xdata[0])*2*np.pi)
    # the problem is there, where I transform the x values.
    xf = np.fft.ifftshift(xf)
    return xf, yf 


# the generated data
x, y = generator(1, 3, 2, delay = 1500, resolution = 0.1) 

# plt.plot(x,y)

xx, yy = ifft_method(x,y)

#if the x values are correctly scaled, the two symmetrical spikes should appear exactly at delay value
# plt.plot(xx, np.abs(yy))

#do the cutting there, which is also irrelevant now

# the problem is there, in fft_method. The x values are not the same as before transforms.
xxx, yyy = fft_method(xx, yy) 
plt.plot(xxx, np.abs(yyy))

#and it should look like this:
#xs = np.linspace(x[0], x[-1], len(x))
#plt.plot(xs, np.abs(yyy))

plt.grid()
plt.show()

A little physical background, if you need explanation:

As in the generator function in the code, I stands for intensity and omega stands for angular frequency. Let's take the Inverse Fourier Transfrom of the generated data:

$F\{I(\omega)\} = F\{I_s(\omega)\} + F\{I_r(\omega)\} + F\{2\cdot\sqrt{I_s(\omega) I_r(\omega)}\cdot \cos({\Phi(\omega)})\}$

where F denotes the Inverse Fourier Transfrom.

Thus, in time domain: $I(t) = I_s(t) + I_r(t) + I_i(t-\tau) + I_i(t+\tau)$

where $I_i(t-\tau)$ and $I_i(t+\tau)$ are called interference terms, and $I_s(t) + I_r(t)$ doesn't matter to us, because they are slowly changing with regards to $\omega$, they appear at 0. Two symmerical spikes should appear at which is the given $\tau$ currently is (called delay in the code above).

My main problem is that I can't figure out how to come back to the original ang. frequency domain.

$\endgroup$
  • $\begingroup$ Since you know the time step, why do you need to extract it from xx? np.arange(start, stop, step), or np.arange(length(xx))*d $\endgroup$ – Dan Boschen Oct 25 at 22:14
  • $\begingroup$ I added more details and reproducible example. $\endgroup$ – Péter Leéh Oct 27 at 11:13
  • $\begingroup$ What about my previous comment? You know the sample rate, you know the length of the data, so you can easily create xf yourself using techniques I described rather than trying to translate fftfreq. $\endgroup$ – Dan Boschen Oct 27 at 11:22
  • $\begingroup$ Oh I misunderstood it first, I'll definitely try that and update the post. $\endgroup$ – Péter Leéh Oct 27 at 11:26

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.