应该如何重新调整 fft 点以获得与解析解相同的结果?
how fft points should be rescaled to get the same results as the analytical solution?
我想使用 numpy fft 包进行快速傅里叶变换,然后我尝试比较解析解和快速傅里叶变换之间的结果,虽然我可以从图中看到我已经做到了曲线相似,很明显比例不同。
我尝试了几个不同版本的频率(angular频率,频率和波数),但是我所有的尝试都没有用,而且在numpy文档中,并不清楚快速傅里叶是如何transform 被精确定义。例如,我想将时间指数函数傅立叶变换到angular频域,f(t)=Exp(-a|t|), F(w)=a/pi*(a²+w²)(根据我们考虑的频率 space,此解析解有多个版本)
def e(t):
return np.exp(-0.5*abs(t))
def F(w):
return 0.5/(np.pi)*(1/(((0.5)**2)+((w)**2)))
t=np.linspace(0,100,1000)
w=np.fft.fftfreq(len(t))
plt.plot(w,F(w),'o',label='F(w)')
plt.legend()
plt.show()
fourier=np.fft.fft(e(t))
plt.plot(w,fourier,'o')
plt.show()
我专门针对频率尝试了上述代码的多种不同变体,但我仍然没有达到 fft 和解析解相似的地步。谁能帮帮我吗?
Fourier transform can be applied to integrable functions such as np.exp(-0.5*abs(t))
. But the Discrete Fourier Transform computes the Fourier transform of periodic signals. See https://dsp.stackexchange.com/questions/26884/about-fourier-transform-of-periodic-signal and What FFTW Really Computes.
因此,长度为T的帧的DFT对应于周期化帧的傅立叶变换。由于帧从0开始,周期右侧的傅立叶变换计算指数衰减:
如您所见,函数np.exp(-0.5*abs(t))
的一半未显示。我们更正它并添加两侧指数衰减的周期增加部分。
我使用频率作为参数:
import matplotlib.pyplot as plt
import numpy as np
def e(t):
return np.exp(-0.5*abs(t))
def F(w):
return 0.5/(np.pi)*(1/(((0.5)**2)+((w)**2)))
def Fc(xi):
#ok , that's sourced from https://en.wikipedia.org/wiki/Fourier_transform ... Square-integrable functions, one-dimensional, line 207
return 2*0.5/(((0.5)**2)+(4*(np.pi*xi)**2))
framelength=100.
nbsample=1000
def ep(t):
#the periodized negative part is added at the end of the frame.
return np.maximum(np.exp(-0.5*abs(t)),np.exp(-0.5*abs(t-framelength)))
t=np.linspace(0,framelength,nbsample, endpoint=False)
#plotting the periodized signal, to see what's happening
ein=ep(t)
tp=np.linspace(0,10*framelength,10*nbsample, endpoint=False)
periodized=np.zeros(10*nbsample)
for i in range(10):
for j in range(nbsample):
periodized[i*nbsample+j]=ein[j]
plt.plot(tp,periodized,'k-',label='periodized frame')
plt.legend()
plt.show()
fourier=np.fft.fft(ep(t))/np.size(ep(t))*framelength
#comparing the mean is useful to check the overall scaling
print np.mean(ep(t))*framelength
print fourier[0]
print Fc(0)
#the frenquencies of the DFT of a frame of length T are 1/T, 2/T ... and negative for the second part of the array.
xi=np.fft.fftfreq(len(t), framelength/len(t))
# comparison between analytical Fourier transform and dft.
plt.plot(xi,Fc(xi),'o',label='F(xi)')
plt.plot(xi,np.real(fourier),'k-', lw=3, color='red', label='DTF')
plt.legend()
plt.show()
结果如下:
对于实验性非周期性信号,当帧被周期化时会出现人为的不连续性。它导致 spectral leakage and windows 用于减弱不连续性及其影响。其中一个潜力 windows,名为泊松 window,是双侧指数衰减!
我想使用 numpy fft 包进行快速傅里叶变换,然后我尝试比较解析解和快速傅里叶变换之间的结果,虽然我可以从图中看到我已经做到了曲线相似,很明显比例不同。
我尝试了几个不同版本的频率(angular频率,频率和波数),但是我所有的尝试都没有用,而且在numpy文档中,并不清楚快速傅里叶是如何transform 被精确定义。例如,我想将时间指数函数傅立叶变换到angular频域,f(t)=Exp(-a|t|), F(w)=a/pi*(a²+w²)(根据我们考虑的频率 space,此解析解有多个版本)
def e(t):
return np.exp(-0.5*abs(t))
def F(w):
return 0.5/(np.pi)*(1/(((0.5)**2)+((w)**2)))
t=np.linspace(0,100,1000)
w=np.fft.fftfreq(len(t))
plt.plot(w,F(w),'o',label='F(w)')
plt.legend()
plt.show()
fourier=np.fft.fft(e(t))
plt.plot(w,fourier,'o')
plt.show()
我专门针对频率尝试了上述代码的多种不同变体,但我仍然没有达到 fft 和解析解相似的地步。谁能帮帮我吗?
Fourier transform can be applied to integrable functions such as np.exp(-0.5*abs(t))
. But the Discrete Fourier Transform computes the Fourier transform of periodic signals. See https://dsp.stackexchange.com/questions/26884/about-fourier-transform-of-periodic-signal and What FFTW Really Computes.
因此,长度为T的帧的DFT对应于周期化帧的傅立叶变换。由于帧从0开始,周期右侧的傅立叶变换计算指数衰减:
np.exp(-0.5*abs(t))
的一半未显示。我们更正它并添加两侧指数衰减的周期增加部分。
我使用频率作为参数:
import matplotlib.pyplot as plt
import numpy as np
def e(t):
return np.exp(-0.5*abs(t))
def F(w):
return 0.5/(np.pi)*(1/(((0.5)**2)+((w)**2)))
def Fc(xi):
#ok , that's sourced from https://en.wikipedia.org/wiki/Fourier_transform ... Square-integrable functions, one-dimensional, line 207
return 2*0.5/(((0.5)**2)+(4*(np.pi*xi)**2))
framelength=100.
nbsample=1000
def ep(t):
#the periodized negative part is added at the end of the frame.
return np.maximum(np.exp(-0.5*abs(t)),np.exp(-0.5*abs(t-framelength)))
t=np.linspace(0,framelength,nbsample, endpoint=False)
#plotting the periodized signal, to see what's happening
ein=ep(t)
tp=np.linspace(0,10*framelength,10*nbsample, endpoint=False)
periodized=np.zeros(10*nbsample)
for i in range(10):
for j in range(nbsample):
periodized[i*nbsample+j]=ein[j]
plt.plot(tp,periodized,'k-',label='periodized frame')
plt.legend()
plt.show()
fourier=np.fft.fft(ep(t))/np.size(ep(t))*framelength
#comparing the mean is useful to check the overall scaling
print np.mean(ep(t))*framelength
print fourier[0]
print Fc(0)
#the frenquencies of the DFT of a frame of length T are 1/T, 2/T ... and negative for the second part of the array.
xi=np.fft.fftfreq(len(t), framelength/len(t))
# comparison between analytical Fourier transform and dft.
plt.plot(xi,Fc(xi),'o',label='F(xi)')
plt.plot(xi,np.real(fourier),'k-', lw=3, color='red', label='DTF')
plt.legend()
plt.show()
结果如下:
对于实验性非周期性信号,当帧被周期化时会出现人为的不连续性。它导致 spectral leakage and windows 用于减弱不连续性及其影响。其中一个潜力 windows,名为泊松 window,是双侧指数衰减!