0

我正在尝试更多地了解噪声、功率谱密度 (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：

  1. 计算前噪声的时间分布

  2. 从时间分布计算的自相关函数（我知道 x 轴比例是错误的，但这对其他任何地方的代码都没有贡献

  3. 功率谱密度，虽然不应该是美丽的对称

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

4

1 回答 1

0

First,your function is written incorrect, you take fft from autocorr which is not your initial signal array. Funny thing, because of the nature of rounding error you will get whit noise like psd from it as well. Second, you calculating frequencies axis wrong, as they should be extended to N/(t*2) (which is Nyquist fr.). Instead, since your freq_axis array length is N, you retrieve N elements from you sig array, and because of that you just read negative frequencies with the same psd, which causes you the symmetry (taking log convert you variable to real and all the Fourier coefficients for negative freqs are just conjugates for positive ones). replacing your code with the following one give a perfectly good result:

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

don't forget to import

import matplotlib.mlab as ml
于 2016-01-29T21:49:08.187 回答