Масштабирование по оси абсцисс после IFFT-FFT - PullRequest
Новая
       15

Масштабирование по оси абсцисс после IFFT-FFT

1 голос
/ 25 октября 2019

Подробности см. В редактировании ниже.

У меня есть набор данных, на котором мне нужно выполнить, и IFFT, вырезать значимую его часть (путем умножения на гауссову кривую), затем БПФ обратно. Первый находится в угловой частотной области, поэтому IFFT приводит к временной области. Тогда обратное БПФ должно снова привести к угловой частоте, но я не могу найти решение, как вернуть исходную область. Конечно, это легко для значений y: 

yf = np.fft.ifft(y)
#cut the valueable part there..
np.fft.fft(yf)

Для преобразований значения x я использую np.fft.fftfreq следующим образом: 

# 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))

Однако выполнение 

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

полностью не возвращает мне первоначальные значения x. Это то, что я неправильно понимаю?

Вопрос можно задать обобщенно, например: 

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 from xx? Is it even possible?

Я пытался использовать временную переменную, где я храню исходные значения x,но это не слишком элегантноЯ ищу какую-то обратную сторону fftfreq и вообще возможное лучшее решение этой проблемы. Спасибо.

РЕДАКТИРОВАТЬ : Я предоставлю код в конце. У меня есть набор данных, который имеет угловую частоту по оси х и интенсивность по оси у. Я хочу выполнить IFFT, чтобы перейти во временную область. К сожалению, значения x не расположены равномерно, поэтому перед IFFT необходима (линейная) интерполяция. Затем во временной области преобразование выглядит следующим образом: See this picture

Следующим шагом является вырезание одного из симметричных пиков с помощью гауссовой кривой, затем БПФ обратно в угловую частотную область (то же самое, где мы начали). Моя проблема в том, что, когда я перевожу ось X для IFFT (что я считаю правильным), я не могу вернуться в исходную угловую частотную область. Вот код, который также включает генератор для набора данных. 

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()
...