MATLAB プログラムを Python に変換しようとしています。クロス パワー スペクトル密度関数の設定に問題があり、Matlab と一致する結果が得られません。
MATLAB コードで使用される関数は次のように記述されます。
[Pxy,f] = cpsd(x,y,M,round(M/2),M,fs);
私のコードで利用可能なドキュメントでは、M = 128 (FFT ポイントの数) および fs = 25.0 (サンプリング周波数 [Hz]) を読みました。x と y は、加速度データの行ベクトル 1x751 です。
[pxy,f] = cpsd(x,y,window,noverlap,f,fs)使用された関数には 6 つの引数があるため、これはプログラマが MATLAB ライブラリから呼び出すことを意図した関数であると仮定します。
この関数は、f で指定された周波数でのクロス パワー スペクトル密度推定値を返します。 (f が周波数として定義されていないのに、変数 M がそこに渡され、それが FFT ポイントの数であることは私を悩ませますが、これは間違いではないと仮定しましょう)。
では、 scipy.signal.csdを使ってこの関数を変換したいのですが、2 つの問題があります。
- ウィンドウは MatLab では整数として定義されていますが、scipy の csd では、ウィンドウをタプル、文字列、または array_like オブジェクトとしてのみ許可しています。
- scipy の csd には、特定の周波数でのクロス パワー スペクトル密度推定値を返すことができる引数はありません。
番号 1 では、ウィンドウを次のように定義し
window = hamming(M, sym=False)
ました。MATLAB の csd でウィンドウが整数として渡されるときに指定された既定のハミング ウィンドウであるため、ハミング ウィンドウを選択します ( 「ウィンドウが整数の場合、cpsd は x と y を長さのウィンドウのセグメントと、その長さのハミング ウィンドウを持つ各セグメントのウィンドウ」 ) であり、スペクトル分析を行っていることを考えると対称にしなかったため、周期的なウィンドウを使用することは理にかなっています。
番号2については、解決策がありません。
これは、Python コードで設定した関数です。
M = 128
fs = 25.0
window = hamming(M, sym=False)
noverlap = np.round(M/2)
f, fxx = signal.csd(x,y,fs=fs,window=window,noverlap=noverlap
結果は、Pxy (クロス パワー スペクトル密度) に関しては一致しませんが、周波数に関しては完全です。これらは、matlab 結果の最初の要素です。
1.3590e-05
3.4354e-05
4.5282e-05
6.2549e-05
5.7965e-05
4.9697e-05
5.5413e-05
これは私がPythonから得たものですが:
2.04688576e-06
3.37540142e-05
4.51821900e-05
6.19997501e-05
5.78926181e-05
5.00058106e-05
5.53681683e-05
次のように、matlplotlib (ドキュメントはこちら) で単純なクロス スペクトル密度関数を使用してみました。
fxx, f = mlab.csd(x,y,NFFT=M,Fs=fs,noverlap=noverlap)
さらに一致する結果が得られますが、まだ完全ではありません。
1.18939608e-05
3.45206717e-05
4.56902859e-05
6.45083475e-05
5.73594952e-05
5.01539145e-05
5.34534933e-05
目的は、変換で発生する可能性のある数値エラーを取り除くことではなく、一致する入力でクロス パワー スペクトル密度を操作することです。
誰でも助けることができますか？よろしくお願いします!!!