Parseval定理不适用于ifft [英] Parseval's theorem doesn't work with ifft

查看:66
本文介绍了Parseval定理不适用于ifft的处理方法,对大家解决问题具有一定的参考价值,需要的朋友们下面随着小编来一起学习吧!

问题描述

我有一个时间信号,然后计算其傅立叶变换以获得频率信号.根据Parseval定理,两个信号具有相同的能量.我成功地用Python演示了它.但是,当我计算频率信号的傅立叶逆变换时,能量不再守恒.这是我的代码:

I have a temporal signal and I calculate its Fourier Transform to get the frequencial signal. According to Parseval's theorem, the two signals have the same energy. I successfully demonstrate it with Python. However, when I calculate the inverse Fourier Transform of the frequencial signal, the energy is no longer conserved. Here is my code:

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))

执行此代码后,您可以看到两个第一个能量是相同的,而第三个能量要比第一个第一个能量小几个数量级,尽管我应该找到相同的能量.这里发生了什么事?

After executing this code, you can see that the two first energies are the same, whereas the third is orders magnitude less than the two first, although I am supposed to find the same energy. What happened here?

推荐答案

与其他软件相比,numpy FFT过程实际上会调整序列长度,这样您就可以得到

The numpy FFT procedures actually and in contrast to other software do adjust for the sequence length, so that you get

nf.ifft(nf.fft(gx)) == gx

直到出现一些浮点错误.如果dxdk是按常规方法计算的,则dk*dx=(2*pi)/N仅适用于未经调整的FFT例程.

up to some floating point error. If your dx and dk are computed the usual way, then dk*dx=(2*pi)/N which only works for unadjusted FFT routines.

您可以使用

测试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可以由

which tells us that the fft is the usual unadjusted transform and ifft divides the result by sequence length N=num. The typical ifft can be emulated by

gxx = (nf.fft(gk.conj())).conj()

那你明白了

gx == gxx/1000

最多浮点错误.或者,您可以使用

up to floating point errors. Or you can reverse the adjustment using

#inverse FT
gx_ = nf.ifft(nf.ifftshift(gkk))*(num*dk)/(2 * np.pi)**0.5

这篇关于Parseval定理不适用于ifft的文章就介绍到这了,希望我们推荐的答案对大家有所帮助,也希望大家多多支持IT屋!

查看全文
登录 关闭
扫码关注1秒登录
发送“验证码”获取 | 15天全站免登陆