Подтвердить что ты не робот

Парасевальная теорема в Python

Я пытаюсь получить некоторое влияние на функциональность Python fft, и одна из странных вещей, на которые я наткнулся, заключается в том, что Парасевальная теорема, похоже, не применяется, так как теперь он имеет значение около 50, тогда как оно должно быть 0.

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

pi = np.pi

tdata = np.arange(5999.)/300
dt = tdata[1]-tdata[0]

datay = np.sin(pi*tdata)+2*np.sin(pi*2*tdata)
N = len(datay)

fouriery = abs(fftpack.rfft(datay))/N

freqs = fftpack.rfftfreq(len(datay), d=(tdata[1]-tdata[0]))

df = freqs[1] - freqs[0]

parceval = sum(datay**2)*dt - sum(fouriery**2)*df
print parceval

plt.plot(freqs, fouriery, 'b-')
plt.xlim(0,3)
plt.show()

Я уверен, что это нормализация, но я, похоже, не могу ее найти, так как вся информация, которую я могу найти об этой функции, - это документация scipy.fftpack.rfft.

Ответ 1

Ваш коэффициент нормализации исходит из попытки применить теорему Парсеваля для преобразования Фурье непрерывного сигнала к дискретной последовательности. На боковой панели статьи в Википедии о дискретном преобразовании Фурье обсуждается связь преобразования Фурье, ряда Фурье, дискретного преобразования Фурье и выборки с гребни Дирака.

Короче говоря, теорема Парсева, применяемая к DFT, не требует интеграции, но суммирование: a 2*pi, которую вы создаете путем многократного на dt и df ваши суммы.

Заметьте также, что, поскольку вы используете scipy.fftpack.rfft, то, что вы получаете, не является надлежащим образом ДПФ ваших данных, но только положительная его половина, так как отрицание будет симметрично ему. Поэтому, поскольку вы добавляете только половину данных, плюс 0 в терминах DC, там отсутствует 2, чтобы добраться до 4*pi, найденного @unutbu.

В любом случае, если datay содержит вашу последовательность, вы можете проверить теорему Парсеваля следующим образом:

fouriery = fftpack.rfft(datay)
N = len(datay)
parseval_1 = np.sum(datay**2)
parseval_2 = (fouriery[0]**2 + 2 * np.sum(fouriery[1:]**2)) / N
print parseval_1 - parseval_2

Используя scipy.fftpack.fft или numpy.fft.fft, второе суммирование не нужно брать на себя странная форма:

fouriery_1 = fftpack.fft(datay)
fouriery_2 = np.fft.fft(datay)
N = len(datay)
parseval_1 = np.sum(datay**2)
parseval_2_1 = np.sum(np.abs(fouriery_1)**2) / N
parseval_2_2 = np.sum(np.abs(fouriery_2)**2) / N
print parseval_1 - parseval_2_1
print parseval_1 - parseval_2_2