Main Content

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)')

Figure contains an axes object. The axes object with title Periodogram Using FFT contains an object of type line.

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

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

Figure contains an axes object. The axes object with title Periodogram Power Spectral Density Estimate contains an object of type line.

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

正規化された周波数を使用した入力

正規化周波数を使用した入力について、fft を使用してピリオドグラムを作成します。正弦波から成り、N(0,1) 加法性ノイズをもつ信号を作成します。この正弦波は π/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)')

Figure contains an axes object. The axes object with title Periodogram Using FFT contains an object of type line.

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

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

Figure contains an axes object. The axes object with title Periodogram Power Spectral Density Estimate contains an object of type line.

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

正規化された周波数を使用した複素数値の入力

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

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

fft を使用してピリオドグラムを求めます。入力が複素数値であるため、[0,2π) ラジアン/サンプルでピリオドグラムを求めます。結果をプロットします。

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)')

Figure contains an axes object. The axes object with title Periodogram Using FFT contains an object of type line.

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

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

Figure contains an axes object. The axes object with title Periodogram Power Spectral Density Estimate contains an object of type line.

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

参考

アプリ

関数