- 周波数帯域の配列インデックスをfind関数を用いて設定しています。
- 周波数帯域のパワースペクトル密度を平均値で算出
- FFTに窓関数を設定していませんが、波形によっては検討が必要な場合もあります。
FFT を使用したパワースペクトル密度推定(PSD)について
124 views (last 30 days)
Show older comments
現在、立っている状態での重心動揺の信号における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
ご教示いただけると幸いです。宜しくお願い致します。
0 Comments
Accepted Answer
covao
on 1 Mar 2024
Edited: covao
on 2 Mar 2024
ご質問のコードを参考に、特定の周波数帯域でのパワースペクトル密度を算出しています。
(生成AIを用いています)
% 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)] );
% 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;
5 Comments
Morio Nishigaki
on 5 Mar 2024
データの切り出し方により、0Hz成分が出力されてしまうことはあり得ます。例えばサイン波で0~540度までの1.5波分を切り出してしまったとすると、プラス側の0Hz成分が多いので、FFTの出力には0Hz成分が出てきます(本来の時間切り出しなしのサイン波では0Hz成分がないにもかからわず)。
簡単に0Hz成分を除去したい場合は、信号Sのとき、S0 = S - mean(S); で求めたS0を使うとよいと思います。
質問の意図から外れていたらすみません。
More Answers (0)
See Also
Categories
Find more on ビッグ データの処理 in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!