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. Power Spectral Denisty、美しく左右対称ではありませんが、

この対称性の原因について誰かが提供できる助けがあれば、最も役に立ちます。単純な正弦波信号でコードをテストしたところ、期待どおりの結果が得られました (対称性なし)。

4

1 に答える 1

0

まず、関数が正しく記述されていません。最初の信号配列ではない autocorr から fft を取得します。面白いことに、丸め誤差の性質により、psd のようなノイズも発生します。次に、周波数軸を N/(t*2) (ナイキスト周波数) に拡張する必要があるため、周波数軸を間違って計算します。代わりに、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 の結果

インポートすることを忘れないでください

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