2013-05-23 3 views
0

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

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) 

ответ

3

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

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

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

myfft[wn:] = 0 

к:

myfft[wn:-wn] = 0 
+0

Отлично - вы совершенно правы. Спасибо, парни. Мне нравится, как быстро люди могут реагировать на SO! – nrob

+0

Обязательно отметьте вопрос как ответ, если вы удовлетворены. – M456

3

Вы удаление половины спектра, когда вы делаете 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 
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) 

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

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

Смежные вопросы