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

Question

我试图通过改变给定数组 (y) 的傅立叶系数的相位来生成一些随机的一维序列 (r)。我假设这两个序列（y 和 r）的振幅谱在执行此操作后不会改变，如果我使用 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。根据文档，这对应于奈奎斯特频率，如果数据大小为奇数，差异不会消失。我不明白是什么导致了差异，或者我应该如何生成具有相同振幅谱的实数值数组。

提前致谢。

Answer 1

我认为 rfft 的（单个）"Nyquist bin" 是问题所在，假设 rfft 使用：

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)

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

Amplitude spectrum change after altering phases of Fourier coefficients

python

random

numpy

fft