FFT を使用したパワースペ​クトル密度推定(PS​D)について

現在、立っている状態での重心動揺の信号におけるPSDを算出するためにFFTの解析をかけております。

指定した周波数帯域のPSDを取得するためにはどうすればよろしいでしょうか？指定する周波数帯域は、0-0.3Hz, 0.3-1.0Hz, 1.0-3.0Hzです。

現在のコードでは、0-0.3Hz帯域でとても大きい値が出てしまっており、原因が分からず質問させていただきました。

下記がコードとなります。

fs = 1000.00;

N = length(time);

fft_COP_X = fft(filt_COP_X);

fft_COP_Y = fft(filt_COP_Y);

fft_COP_X = fft_COP_X(1:N/2+1);

p_fft_COP_X = (1/(fs*N)) * abs(fft_COP_X).^2;

p_fft_COP_X(2:end-1) = 2*p_fft_COP_X(2:end-1);

freq = 0:fs/length(time):fs/2;

fft_COP_Y = fft_COP_Y(1:N/2+1);

p_fft_COP_Y = (1/(fs*N)) * abs(fft_COP_Y).^2;

p_fft_COP_Y(2:end-1) = 2*p_fft_COP_Y(2:end-1);

freq = 0:fs/length(time):fs/2;

Plow_fft_COP_X = sum(p_fft_COP_X(1:10));%0-0.3Hz 低周波数帯域 Power

Pmedium_fft_COP_X = sum(p_fft_COP_X(10:32));%0.3-1.0Hz 中周波数帯域 Power

Phigh_fft_COP_X = sum(p_fft_COP_X(32:92));%1.0-3.0Hz 高周波数帯域 Power

Plow_fft_COP_Y = sum(p_fft_COP_Y(1:10));%0-0.3Hz 低周波数帯域 Power

Pmedium_fft_COP_Y = sum(p_fft_COP_Y(10:32));%0.3-1.0Hz 中周波数帯域 Power

Phigh_fft_COP_Y = sum(p_fft_COP_Y(32:92));%1.0-3.0Hz 高周波数帯域 Power

ご教示いただけると幸いです。宜しくお願い致します。

Answer 1

[covao]

ご質問のコードを参考に、特定の周波数帯域でのパワースペクトル密度を算出しています。

（生成AIを用いています）

• 周波数帯域の配列インデックスをfind関数を用いて設定しています。
• 周波数帯域のパワースペクトル密度を平均値で算出
• FFTに窓関数を設定していませんが、波形によっては検討が必要な場合もあります。

% Define the sampling frequency and time vector
fs = 1000.0;
t = 0:1/fs:20;

% Generate the signals (combination of 0.15Hz, 0.65Hz, and 2Hz)
filt_COP_X = sin(2*pi*0.15*t) + 0.5*sin(2*pi*0.65*t) + 0.2*sin(2*pi*2*t);

% Perform FFT analysis
N = length(filt_COP_X);
fft_COP_X = fft(filt_COP_X);
fft_COP_X = fft_COP_X(1:floor(N/2)+1);
p_fft_COP_X = (1/(fs*N)) * abs(fft_COP_X).^2;
p_fft_COP_X(2:end-1) = 2*p_fft_COP_X(2:end-1);
freq = linspace(0, fs/2, length(p_fft_COP_X));

% Calculate indices corresponding to specified frequency bands
idx_low = find(freq >= 0 & freq <= 0.3);
idx_medium = find(freq > 0.3 & freq <= 1.0);
idx_high = find(freq > 1.0 & freq <= 3.0);

% Calculate power in each frequency band
Plow_fft_COP_X = mean(p_fft_COP_X(idx_low));
Pmedium_fft_COP_X = mean(p_fft_COP_X(idx_medium));
Phigh_fft_COP_X = mean(p_fft_COP_X(idx_high));


% Output the results
disp(['PSD X: Low=', num2str(Plow_fft_COP_X), ', Mid=', num2str(Pmedium_fft_COP_X), ', High=', num2str(Phigh_fft_COP_X)] );

PSD X: Low=1.4286, Mid=0.17856, High=0.0099974

% Plot the power spectrum and original signals
figure;
% Original Signals
subplot(2, 1, 1);
plot(t, filt_COP_X);
title('Original Signal COP X');
xlabel('Time (s)');
ylabel('Amplitude');
grid on;

% Power Spectrum
subplot(2, 1, 2);
plot(freq, p_fft_COP_X);
title('Power Spectrum COP X');
xlabel('Frequency (Hz)');
ylabel('Power');
xlim([0, 3]); % Display up to 3Hz
grid on;

[拓也 矢田]

お忙しいところ、ご対応ありがとうございました。

大変参考になりました。

[拓也 矢田]

平素よりお世話になっております。

教えていただいたコードでPSDを算出したのですが、全てのデータ(異なる試技)で同じ値が算出されてしまいます。MATLAB、FFTともに初学者のため、原因が分かりません。原因が分かる方がいらっしゃれば、ご教示いただけると幸いです。宜しくお願い致します。

[covao]

各データの差が小さいためか、データ処理上の問題か、切り分けが必要と考えられます。

波形と周波数特性をプロットして、明らかに違いがあるようでしたコードに原因があると思われます。

[Morio Nishigaki]

データの切り出し方により、0Hz成分が出力されてしまうことはあり得ます。例えばサイン波で0～540度までの1.5波分を切り出してしまったとすると、プラス側の0Hz成分が多いので、FFTの出力には0Hz成分が出てきます（本来の時間切り出しなしのサイン波では0Hz成分がないにもかからわず）。

簡単に0Hz成分を除去したい場合は、信号Sのとき、S0 = S - mean(S); で求めたS0を使うとよいと思います。

質問の意図から外れていたらすみません。

[拓也 矢田]

お忙しいところ、ご対応いただきありがとうございました。S0 = S - mean(S); で求めたS0を用いることによって、解決しました。大変勉強になりました。ありがとうございました。