17

我试图掌握 Python 的 fft 功能，而我偶然发现的一件奇怪的事情是Parseval 定理似乎并不适用，因为它现在给出了大约 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 文档。

4

1 回答 1

24

您的归一化因子来自尝试将 Parseval 定理应用于连续信号到离散序列的傅里叶变换。在关于离散傅里叶变换的维基百科文章的侧板上，有一些关于傅里叶变换、傅里叶级数、离散傅里叶变换和狄拉克梳采样的关系的讨论。

长话短说，当应用于 DFT 时，Parseval 定理不需要积分，而是求和：a2*pi您通过乘以dt和df求和来创建。

另请注意，因为您使用的是 ，所以您scipy.fftpack.rfft得到的不是数据的正确 DFT，而只是数据的正半部分，因为负部分与它对称。因此，由于您只添加了一半的数据，加上0DC 术语中的数据，因此缺少2找到4*pi@unutbu 找到的数据。

在任何情况下，如果datay持有您的序列，您可以验证 Parseval 定理如下：

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.fftornumpy.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
于 2012-12-23T16:18:38.107 回答