ドキュメンテーション

最新のリリースでは、このページがまだ翻訳されていません。 このページの最新版は英語でご覧になれます。

FFT を使用したパワー スペクトル密度推定

この例では、fft を使用して、ピリオドグラムと等価なノンパラメトリック PSD (パワー スペクトル密度) 推定を求める方法を示します。例においては、偶数長の入力、正規化された周波数、片側および両側 PSD 推定に対し、fft の出力を適切にスケーリングする方法が示されます。fftperiodogram の両方を使用して、1 kHz でサンプリングされた偶数長の信号についてピリオドグラムを求めます。結果を比較します。

100 Hz 正弦波から成り、N(0,1) 加法性ノイズをもつ信号を作成します。サンプリング周波数は 1 kHz です。信号長は 1000 サンプルです。結果に再現性をもたせるために、乱数発生器の既定の設定を使用します。

rng default
Fs = 1000;
t = 0:1/Fs:1-1/Fs;
x = cos(2*pi*100*t) + randn(size(t));

fft を使用してピリオドグラムを求めます。信号は実数値で、偶数長です。信号が実数値なので、正または負の周波数に対するパワー推定だけを必要とします。合計したパワーを保存するために、両セット (正および負の周波数) に現れているすべての周波数に係数 2 を乗算します。周波数 0 (DC) およびナイキスト周波数は、2 回は現れません。結果をプロットします。

N = length(x);
xdft = fft(x);
xdft = xdft(1:N/2+1);
psdx = (1/(Fs*N)) * abs(xdft).^2;
psdx(2:end-1) = 2*psdx(2:end-1);
freq = 0:Fs/length(x):Fs/2;

plot(freq,10*log10(psdx))
grid on
title('Periodogram Using FFT')
xlabel('Frequency (Hz)')
ylabel('Power/Frequency (dB/Hz)')

periodogram を使用してピリオドグラムを計算してプロットします。2 つの結果が同じであることを示します。

periodogram(x,rectwin(length(x)),length(x),Fs)

mxerr = max(psdx'-periodogram(x,rectwin(length(x)),length(x),Fs))
mxerr =

   3.4694e-18

正規化された周波数を使用した入力について、fft を使用してピリオドグラムを作成します。正弦波から成り、N(0,1) 加法性ノイズをもつ信号を作成します。この正弦波は $\pi/4$ ラジアン/サンプルの角周波数をもちます。結果に再現性をもたせるために、乱数発生器の既定の設定を使用します。

rng default
n = 0:999;
x = cos(pi/4*n) + randn(size(n));

fft を使用してピリオドグラムを求めます。信号は実数値で、偶数長です。信号が実数値なので、正または負の周波数に対するパワー推定だけを必要とします。合計したパワーを保存するために、両セット (正および負の周波数) に現れているすべての周波数に係数 2 を乗算します。周波数 0 (DC) およびナイキスト周波数は、2 回は現れません。結果をプロットします。

N = length(x);
xdft = fft(x);
xdft = xdft(1:N/2+1);
psdx = (1/(2*pi*N)) * abs(xdft).^2;
psdx(2:end-1) = 2*psdx(2:end-1);
freq = 0:(2*pi)/N:pi;

plot(freq/pi,10*log10(psdx))
grid on
title('Periodogram Using FFT')
xlabel('Normalized Frequency (\times\pi rad/sample)')
ylabel('Power/Frequency (dB/rad/sample)')

periodogram を使用してピリオドグラムを計算してプロットします。2 つの結果が同じであることを示します。

periodogram(x,rectwin(length(x)),length(x))

mxerr = max(psdx'-periodogram(x,rectwin(length(x)),length(x)))
mxerr =

   1.4211e-14

正規化された周波数をもつ複素数値の入力について、fft を使用してピリオドグラムを作成します。信号は、複素数値の N(0,1) ノイズを伴う $\pi/4$ ラジアン/サンプルの角周波数をもつ複素指数です。再現性のある結果を得るために、乱数発生器を既定の状態に設定します。

rng default
n = 0:999;
x = exp(1j*pi/4*n) + [1 1j]*randn(2,length(n))/sqrt(2);

fft を使用してピリオドグラムを求めます。入力が複素数値なので、 $(-\pi,\pi]$ ラジアン/サンプルでピリオドグラムを求めます。結果をプロットします。

N = length(x);
xdft = fft(x);
psdx = (1/(2*pi*N)) * abs(xdft).^2;
freq = 0:(2*pi)/N:2*pi-(2*pi)/N;

plot(freq/pi,10*log10(psdx))
grid on
title('Periodogram Using FFT')
xlabel('Normalized Frequency (\times\pi rad/sample)')
ylabel('Power/Frequency (dB/rad/sample)')

periodogram を使用してピリオドグラムを求めてプロットします。PSD 推定を比較します。

periodogram(x,rectwin(length(x)),length(x),'twosided')

mxerr = max(psdx'-periodogram(x,rectwin(length(x)),length(x),'twosided'))
mxerr =

   2.8422e-14

この情報は役に立ちましたか?