改变傅里叶系数相位后振幅谱的变化

2024-03-29 13:00:45 发布

您现在位置:Python中文网/ 问答频道 /正文

我试图通过改变给定数组的Fourier系数的相位来生成一些随机的1d序列(r)。我假设这两个序列(yr)的振幅谱在这样做之后不会改变,如果我使用numpy.fft.fft(),它们确实不会改变,就像这样:

import numpy as np
import numpy.fft as fft

n=512
x=np.linspace(0,10*np.pi,n)
y=np.sin(x)+2*np.sin(2*x)+2*np.sin(5*x)

#-------------Get Fourier coefficients-------------
cf1=fft.fft(y)
amp1=np.abs(cf1)
theta1=np.angle(cf1)

#-Randomly alter phase keeping amplitude unchanged-
theta2=np.random.normal(0,1.,size=theta1.shape)+theta1
cf2=amp1*np.cos(theta2)+1j*amp1*np.sin(theta2)

#-----------------------IFFT-----------------------
y2=fft.ifft(cf2)

#------------Compare amplitude spectrum------------
cf3=fft.fft(y2)
amp2=np.abs(cf3)

import matplotlib.pyplot as plt
figure=plt.figure()
ax=figure.add_subplot(111)
ax.plot(amp1-amp2,'k-')

plt.show(block=False)

结果图是1e-13顺序的随机序列。如此有效地保持不变。在

但我对生成随机的真实数据而不是复杂的数据感兴趣。使用numpy.fft.rfft()numpy.fft.irfft(),除了最后一个频率(amp1[-1]amp2[-1])之外,所有频率的振幅都一致,并且差值为0.1级。根据文档,这与Nyquist频率相对应,如果数据大小为奇数,则差异不会消失。我不知道是什么造成了这种差异,也不知道如何生成一个具有相同振幅谱的实值阵列。在

提前谢谢。在


Tags: importfftnumpyasnpplt序列sin
1条回答
网友
1楼 · 发布于 2024-03-29 13:00:45

我认为rfft的(单个)“Nyquist-bin”就是问题所在,在{}使用的假设下,您不应该改变它的相位:

When the input is purely real, its transform is Hermitian, i.e., the component at frequency f_k is the complex conjugate of the component at frequency -f_k, which means that for real inputs there is no information in the negative frequency components that is not already available from the positive frequency components. The family of rfft functions is designed to operate on real inputs, and exploits this symmetry by computing only the positive frequency components, up to and including the Nyquist frequency. Thus, n input points produce n/2+1 complex output points. The inverses of this family assumes the same symmetry of its input, and for an output of n points uses n/2+1 input points.

当我修补最后一个rfft bin以保持原始相位时,图看起来不错

cf1=fft.rfft(y)
amp1=np.abs(cf1)
theta1=np.angle(cf1)
tlast=theta1[-1] # quick patch, save last (Nyquist) phase
#-Randomly alter phase keeping amplitude unchanged-
theta2=np.random.normal(0,1.,size=theta1.shape)+theta1
theta2[-1] = tlast # restore Nyquist bin original phase
cf2=amp1*np.cos(theta2)+1j*amp1*np.sin(theta2)

相关问题 更多 >