ご質問のコードを参考に、特定の周波数帯域でのパワースペクトル密度を算出しています。
（生成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)] );
% 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;
