问题 Parseval在Python中的定理


我试图掌握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文档


6785
2017-12-23 13:54


起源



答案:


你的 归一化因子 试图将Parseval定理应用于连续信号的傅立叶变换到离散序列。在侧面板上 关于离散傅立叶变换的维基百科文章 对傅里叶变换,傅立叶级数,离散傅里叶变换和采样的关系进行了讨论 狄拉克梳理

使长话短说, Parseval定理,当应用于DFT时,不需要整合,而是总结:a 2*pi 你是通过乘法来创造的 dt 和 df 你的总结。

另请注意,因为您正在使用 scipy.fftpack.rfft,你得到的不是数据的DFT,而只是正数的一半,因为负数将与它对称。所以,因为你只添加了一半的数据,加上 0在DC术语中,缺少了 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.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

15
2017-12-23 16:18



谢谢你的纠正! - unutbu
要注意:在某种程度上,这是一般问题的一个方面,浮点数与实数不同。 - Marcin
@Marcin是的,改了 == 到了 - 在检查... - Jaime
@Jamie - 很好的解释!另外,您可以使用 np.allclose 检查浮动的大致相等而不是打印差异。例如。 assert np.allclose(parseval_1, parseval_2_1) (allclose 主要用于检查两个数组中的所有项是否相似,因此名称,但它适用于标量。) - Joe Kington


答案:


你的 归一化因子 试图将Parseval定理应用于连续信号的傅立叶变换到离散序列。在侧面板上 关于离散傅立叶变换的维基百科文章 对傅里叶变换,傅立叶级数,离散傅里叶变换和采样的关系进行了讨论 狄拉克梳理

使长话短说, Parseval定理,当应用于DFT时,不需要整合,而是总结:a 2*pi 你是通过乘法来创造的 dt 和 df 你的总结。

另请注意,因为您正在使用 scipy.fftpack.rfft,你得到的不是数据的DFT,而只是正数的一半,因为负数将与它对称。所以,因为你只添加了一半的数据,加上 0在DC术语中,缺少了 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.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

15
2017-12-23 16:18



谢谢你的纠正! - unutbu
要注意:在某种程度上,这是一般问题的一个方面,浮点数与实数不同。 - Marcin
@Marcin是的,改了 == 到了 - 在检查... - Jaime
@Jamie - 很好的解释!另外,您可以使用 np.allclose 检查浮动的大致相等而不是打印差异。例如。 assert np.allclose(parseval_1, parseval_2_1) (allclose 主要用于检查两个数组中的所有项是否相似,因此名称,但它适用于标量。) - Joe Kington