Восстановить оригинальный сигнал с БПФ в python - PullRequest
Новая
       21

Восстановить оригинальный сигнал с БПФ в python

4 голосов
/ 13 февраля 2020

У меня есть аналитически сгенерированный спектр, где ось x представляет angular частоту, y представляет интенсивность. Спектр сосредоточен вокруг некоторого значения частоты, которое часто называют центральной частотой сигнала (синий график на рисунке). Я хочу выполнить IFFT для данных во временной области, вырезать их полезную часть с помощью гауссовой кривой, а затем FFT вернуться к исходной области. Моя проблема в том, что после IFFT (FFT (сигнал)) центральная частота теряется: я возвращаю спектр по форме, но он всегда центрируется около 0 (оранжевый график). Spectrum В настоящее время мое решение довольно плохое: я кеширую исходную ось x и восстанавливаю ее после вызовов FFT. Это, очевидно, имеет много недостатков, и я хочу улучшить его. Ниже я включил небольшую демонстрацию, которая демонстрирует проблему. Мой вопрос: можно ли решить это более элегантным способом? Есть ли способ, чтобы центральная частота не терялась во время процесса? 

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

1 Ответ

2 голосов
/ 13 февраля 2020

Я думаю, что fftfreq не делает то, что вы думаете, что делает. xf для fft(ifft(y) идентично x, вам не следует пытаться пересчитать его. Ось X не изменяется при переходе в другую область и затем обратно.

Также обратите внимание, что fftfreq возвращает координаты в частотной области для дискретного преобразования Фурье сигнала данного длина и с заданным интервалом выборки. Он не делает обратное, вы не можете использовать его для определения координат в пространственной области после применения обратного дискретного преобразования Фурье. (Интервал, который он возвращает, действителен, но набор координат - нет.) 

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

Другая проблема с вашим кодом состоит в том, что ifft(y) предполагает фиксированный набор значений вдоль оси x , Ваш x не соответствует этому. Таким образом, сигнал пространственной области, который вы получаете, не имеет смысла.

Запустив ваш код, я вижу, что x работает с 3,0 до 1,0 с шагом 0,0004777. Вам нужно будет дополнить ваши данные так, чтобы значения варьировались от 0,0 до 6,0, при этом область (3.0, 6.0) является сопряженной симметричной c копией региона (0.0, 3.0). Эта область соответствует отрицательным частотам в соответствии с периодичностью частотной области (F [n] == F [n + N], где N - количество выборок). Заполните область (0.0, 1.0) нулями.

Учитывая эту стандартизированную ось x в частотной области, xf = fftfreq(len(xt), d=(xt[1]-xt[0])) должен восстановить ось x. Но вам нужно вычислить xt соответственно: xt = np.linspace(0, 1/(x[1]-x[0]), len(x), endpoint=False) (с x стандартизированной осью частоты DFT).

...