Parseval 定理不适用于 ifft
Parseval's theorem doesn't work with ifft
我有一个时间信号,我计算它的傅立叶变换以获得频率信号。根据帕塞瓦尔定理,这两个信号具有相同的能量。我用 Python 成功地演示了它。但是,当我计算频率信号的傅里叶逆变换时,能量不再守恒。这是我的代码:
import numpy as np
import numpy.fft as nf
import matplotlib.pyplot as plt
#create a gaussian as a temporal signal
x = np.linspace(-10.0,10.0,num=1000)
dx = x[1]-x[0]
sigma = 0.4
gx = (1.0/(2.0*np.pi*sigma**2.0)**0.5)*np.exp(-0.5*(x/sigma)**2.0)
#calculate the spacing of the frequencial signal
f=nf.fftshift(nf.fftfreq(1000,dx))
kk = f*(2.0*np.pi)
dk = kk[1]-kk[0]
#calculate the frequencial signal (FT)
#the convention used here allows to find the same energy
gkk = nf.fftshift(nf.fft(nf.fftshift(gx)))*(dx/(2.0*np.pi)**0.5)
#inverse FT
gx_ = nf.ifft(nf.ifftshift(gkk))*dk/(2 * np.pi)**0.5
#Parseval's theorem
print("Total energy in time domain = "+str(sum(abs(gx)**2.0)*dx))
print("Total energy in freq domain = "+str(sum(abs(gkk)**2.0)*dk))
print("Total energy after iFT = "+str(sum(abs(gx_)**2.0)*dx))
执行这段代码后,可以看到前两个能量相同,而第三个比前两个小一个数量级,虽然我应该找到相同的能量。这里发生了什么?
numpy
FFT 程序实际上与其他软件相比会针对序列长度进行调整,因此您可以得到
nf.ifft(nf.fft(gx)) == gx
最多一些浮点错误。如果您的 dx
和 dk
是按照通常的方式计算的,那么 dk*dx=(2*pi)/N
仅适用于未经调整的 FFT 例程。
您可以使用
测试numpy.fft
的行为
In [20]: sum(abs(gx)**2.0)
Out[20]: 35.226587122763036
In [21]: gk = nf.fft(gx)
In [22]: sum(abs(gk)**2.0)
Out[22]: 35226.587122763049
In [23]: sum(abs(nf.ifft(gk))**2.0)
Out[23]: 35.226587122763014
告诉我们 fft
是通常的未调整变换,ifft
将结果除以序列长度 N=num
。典型的 ifft
可以被
模拟
gxx = (nf.fft(gk.conj())).conj()
然后你就明白了
gx == gxx/1000
多达浮点错误。或者您可以使用
反转调整
#inverse FT
gx_ = nf.ifft(nf.ifftshift(gkk))*(num*dk)/(2 * np.pi)**0.5
我有一个时间信号,我计算它的傅立叶变换以获得频率信号。根据帕塞瓦尔定理,这两个信号具有相同的能量。我用 Python 成功地演示了它。但是,当我计算频率信号的傅里叶逆变换时,能量不再守恒。这是我的代码:
import numpy as np
import numpy.fft as nf
import matplotlib.pyplot as plt
#create a gaussian as a temporal signal
x = np.linspace(-10.0,10.0,num=1000)
dx = x[1]-x[0]
sigma = 0.4
gx = (1.0/(2.0*np.pi*sigma**2.0)**0.5)*np.exp(-0.5*(x/sigma)**2.0)
#calculate the spacing of the frequencial signal
f=nf.fftshift(nf.fftfreq(1000,dx))
kk = f*(2.0*np.pi)
dk = kk[1]-kk[0]
#calculate the frequencial signal (FT)
#the convention used here allows to find the same energy
gkk = nf.fftshift(nf.fft(nf.fftshift(gx)))*(dx/(2.0*np.pi)**0.5)
#inverse FT
gx_ = nf.ifft(nf.ifftshift(gkk))*dk/(2 * np.pi)**0.5
#Parseval's theorem
print("Total energy in time domain = "+str(sum(abs(gx)**2.0)*dx))
print("Total energy in freq domain = "+str(sum(abs(gkk)**2.0)*dk))
print("Total energy after iFT = "+str(sum(abs(gx_)**2.0)*dx))
执行这段代码后,可以看到前两个能量相同,而第三个比前两个小一个数量级,虽然我应该找到相同的能量。这里发生了什么?
numpy
FFT 程序实际上与其他软件相比会针对序列长度进行调整,因此您可以得到
nf.ifft(nf.fft(gx)) == gx
最多一些浮点错误。如果您的 dx
和 dk
是按照通常的方式计算的,那么 dk*dx=(2*pi)/N
仅适用于未经调整的 FFT 例程。
您可以使用
测试numpy.fft
的行为
In [20]: sum(abs(gx)**2.0)
Out[20]: 35.226587122763036
In [21]: gk = nf.fft(gx)
In [22]: sum(abs(gk)**2.0)
Out[22]: 35226.587122763049
In [23]: sum(abs(nf.ifft(gk))**2.0)
Out[23]: 35.226587122763014
告诉我们 fft
是通常的未调整变换,ifft
将结果除以序列长度 N=num
。典型的 ifft
可以被
gxx = (nf.fft(gk.conj())).conj()
然后你就明白了
gx == gxx/1000
多达浮点错误。或者您可以使用
反转调整#inverse FT
gx_ = nf.ifft(nf.ifftshift(gkk))*(num*dk)/(2 * np.pi)**0.5