обратная функция БПФ не совпадает с исходной функцией

Я не понимаю, почему ifft(fft(myFunction)) не совпадает с моей функцией. Кажется, это та же форма, но с коэффициентом 2 (без учета постоянного смещения по оси y). Во всей документации, которую я вижу, говорится, что существует некоторая нормализация, которую не выполняет fft, но ifft должен об этом позаботиться. Ниже приведен пример кода — вы можете видеть, где я использовал коэффициент 2, чтобы дать мне правильный ответ. Спасибо за любую помощь - это сводит меня с ума.

import numpy as np
import scipy.fftpack as fftp
import matplotlib.pyplot as plt
import matplotlib.pyplot as plt

def fourier_series(x, y, wn, n=None):
    # get FFT
    myfft = fftp.fft(y, n)
    # kill higher freqs above wavenumber wn
    myfft[wn:] = 0
    # make new series
    y2 = fftp.ifft(myfft).real
    # find constant y offset
    myfft[1:]=0
    c = fftp.ifft(myfft)[0]
    # remove c, apply factor of 2 and re apply c
    y2 = (y2-c)*2 + c

    plt.figure(num=None)
    plt.plot(x, y, x, y2)
    plt.show()

if __name__=='__main__':
    x = np.array([float(i) for i in range(0,360)])
    y = np.sin(2*np.pi/360*x) + np.sin(2*2*np.pi/360*x) + 5

    fourier_series(x, y, 3, 360)

numpy scipy fft ifft
person nrob    schedule 23.05.2013    source источник

Ответы (2)


arrow_upward
5
arrow_downward

Вы убиваете отрицательные частоты между 0 и -wn.

Я думаю, что вы хотите установить myfft на 0 для всех частот за пределами [-wn, wn].

Измените следующую строку:

myfft[wn:] = 0

to:

myfft[wn:-wn] = 0
person M456    schedule 23.05.2013
comment
Отлично - вы совершенно правы. Спасибо, парни. Мне нравится, как быстро люди могут отвечать на SO! - person nrob; 23.05.2013
comment
Не забудьте отметить вопрос как ответ, если вы удовлетворены. - person M456; 23.05.2013

arrow_upward
8
arrow_downward

Вы удаляете половину спектра, когда делаете myfft[wn:] = 0. Отрицательные частоты находятся в верхней половине массива и являются обязательными.

У вас есть вторая выдумка, чтобы получить ваши результаты, которые берут действительную часть, чтобы найти y2: y2 = fftp.ifft(myfft).real (fftp.ifft(myfft) имеет незначительную мнимую часть из-за асимметрии в спектре).

Исправьте это с помощью myfft[wn:-wn] = 0 вместо myfft[wn:] = 0 и удалите выдумки. Таким образом, фиксированный код выглядит примерно так:

import numpy as np
import scipy.fftpack as fftp
import matplotlib.pyplot as plt    

def fourier_series(x, y, wn, n=None):
    # get FFT
    myfft = fftp.fft(y, n)
    # kill higher freqs above wavenumber wn
    myfft[wn:-wn] = 0
    # make new series
    y2 = fftp.ifft(myfft)

    plt.figure(num=None)
    plt.plot(x, y, x, y2)
    plt.show()

if __name__=='__main__':
    x = np.array([float(i) for i in range(0,360)])
    y = np.sin(2*np.pi/360*x) + np.sin(2*2*np.pi/360*x) + 5

    fourier_series(x, y, 3, 360)

Действительно стоит обратить внимание на промежуточные массивы, которые вы создаете при попытке выполнить обработку сигнала. Неизменно есть подсказки относительно того, что идет не так, что должно направить вас к проблеме. В этом случае вы, взяв на себя настоящую роль, замаскировали проблему и усложнили себе задачу.

Просто добавим еще один быстрый момент: иногда взять действительную часть результирующего массива — это именно то, что нужно сделать. Часто бывает так, что вы получаете мнимую часть выходного сигнала, которая просто связана с числовыми ошибками на входе обратного БПФ. Обычно это проявляется в виде очень маленьких мнимых значений, поэтому действительная часть — это в основном тот же массив.

person Henry Gomersall    schedule 23.05.2013