计算白噪声的功率谱密度时无法解释的对称性答案

【问题标题】：Unexplained symmetry when computing Power Spectral Density of white noise计算白噪声的功率谱密度时无法解释的对称性
【发布时间】：2014-06-19 11:09:43
【问题描述】：

我正在尝试了解有关噪声、功率谱密度 (PSD) 和统计方差的更多信息。关于这一点，我正在尝试计算白噪声的功率谱密度，但是，当我这样做时，我得到了一个非常奇怪的对称性。我的频谱似乎围绕中心频率值对称，这显然是不正确的。我是使用自相关和功率谱密度的新手，所以如果有人能帮助我朝着错误的方向前进，我将不胜感激。

计算PSD的代码：

import numpy as np 
from math import sin, pi, log10
from allan_noise import white,pink,brown, violet
import acor
import numpy as np


#METHOD ONE: AUTOCORRELATION FUNCTION METHOD
def psd1(time_dats):
    #auto_corr = acor.function(time_dats)
    auto_corr = np.correlate(time_dats,time_dats,mode="same")

    #To check autocorrelation
    for i in range(len(auto_corr)):
        print i*.0000001 ,auto_corr[i]

    return fft(auto_corr) 

#DEFINE VARIABLES
t = .0001       #simulation time
N = 100000  #number of data points 
dt = t/N    #time interval between data points

#CREATE SIGNAL
sig = white(N)
df = 1.0/(t)
freq_axis = np.arange(0,N/t,df)

spec = psd1(sig)

#OPEN UP OUTPUT FILE
f = open('data/psdtest_f','w')
g = open('data/psdtest_t','w')

#PRINT OUT DATA
for i in range(N):
    f.write('%g %g\n' %(freq_axis[i],log10(abs(spec[i]))))
    g.write('%g %g\n' %(i*dt, sig[i]))

使用此代码，我生成了以下图表，可以在此处访问https://drive.google.com/#folders/0B8BdiIhqTYw7Vk1EaDFMQW84RHM：

计算前噪声的时间分布
根据时间曲线计算的自相关函数（我知道 x 轴刻度是错误的，但这对其他任何地方的代码都没有影响
功率谱密度，虽然不应该是漂亮的对称

任何人都可以就导致这种对称性的原因提供任何帮助将是最有帮助的。我用一个简单的正弦波信号测试了代码，我得到了预期的结果（没有对称性）

【问题讨论】：

这个问题似乎是题外话，因为它是关于结果的解释，而不是代码实现
它可能更适合例如Signal Processing 或者 Cross-Validated
道歉，虽然我觉得它的代码实现我错了。我不完全理解 numpy 函数是如何工作的。不过如果你觉得它跑题了，我当然可以移动它
是不是你觉得实现不正确？你期待的结果是什么？您是否测试过代码 - 哪些有效，哪些无效？解释代码的请求在这里是题外话，一般 “出了什么问题？” 很难回答 - 请参阅 How to Ask 并考虑提供 minimal example
你说得对，我是 stackoverflow 的新手，而且还在习惯它。我不是要求任何人解释 numpy 库是如何工作的，我只是简单地说我认为这可能是我的错误所在，并且有更多经验的人可能会接受它。我试图在底部说明我已经使用 sin 函数测试了代码并产生了对应于单个频率的峰值。白噪声应该给出平坦的功率谱密度。基本上它应该看起来像我提供的图表的前半部分，让我相信我在某处镜像数据。

标签： python signals symmetric

【解决方案1】：

首先，您的函数编写不正确，您从 autocorr 中获取 fft，这不是您的初始信号数组。有趣的是，由于舍入误差的性质，你也会从中得到像 psd 这样的噪音。其次，您计算频率轴错误，因为它们应该扩展到 N/(t*2) （这是 Nyquist fr.）。相反，由于您的 freq_axis 数组长度为 N，您从 sig 数组中检索 N 个元素，因此您只需使用相同的 psd 读取负频率，这会导致对称性（将 log 转换为实数和所有傅立叶负频率的系数只是正频率的共轭）。用以下代码替换您的代码会产生非常好的结果：

sig = white(N,1,N/t)
(siglog,freq_axis)=ml.psd(sig,N,(N/t), detrend=ml.detrend_none,
   window=ml.window_hanning, noverlap=0, pad_to=None,
    sides='default', scale_by_freq=None)
plt.plot(freq_axis,np.log10(siglog))
plt.show()

matplotlib.mlab.psd result

别忘了导入

import matplotlib.mlab as ml

【讨论】：