Fix the code calculate PSD

1 次查看（过去 30 天）
Tran
Tran 2022-12-16
回答： Star Strider 2022-12-16
I run this code but don't know why the output is very large (while PSD should be very small). Can anyone help me fix this code?
%load data
load('data.mat')
EEG_data = data_par1(:,1);
%design low/high pass
[all,f1] = fir1(9,[0.0625 0.390625],"bandpass",hamming(10));% 4-25Hz
[th,f2] = fir1(9,[0.0625 0.109375],"bandpass",hamming(10)); % 4-7Hz
[a,f3] = fir1(9,[0.109375 0.1875],"bandpass",hamming(10)); % 7-12hz
[b,f4] = fir1(9,[0.1875 0.390625],"bandpass",hamming(10)); % 12 -25Hz
% filtering the data
All = filtfilt(all,f1,EEG_data);
fTheta=filtfilt(th,f2,All);
fAlpha=filtfilt(a,f3,All);
fBeta=filtfilt(b,f4,All);
% calculate PSD
fs =128;
Theta = cpsd(fTheta,fTheta,5*fs,2.5*fs,5*fs,fs);
Alpha = cpsd(fAlpha,fAlpha,5*fs,2.5*fs,5*fs,fs);
Beta = cpsd(fBeta,fBeta,5*fs,2.5*fs,5*fs,fs);
% calculate the average PSD of each signal
psdth = mean(Theta);
psda = mean(Alpha);
psdb = mean(Beta);
  1 个评论
Tran
Tran 2022-12-16
data_par1(:,1) mean that i choose the 1 channels(AF3) to do

请先登录，再进行评论。

请先登录，再回答此问题。

回答（1 个）

Star Strider
Star Strider 2022-12-16
Ran in:
I am not certain what you are doing. The third argument to cpsd is supposed to be a window function, however in your calls to it, it is a scalar, ‘5*fs’.
Also, an order of 9 is quite short for an FIR ffilter. Usually, orders of greater than 30 are required for decent results. Use the freqz function to examine the filter Bode plots.
I use elliptic filters here, courtesy of the bandpass function (although designing them from command-line function calls is straightforward) —
LD = load(websave('data','https://www.mathworks.com/matlabcentral/answers/uploaded_files/1233172/data.mat'))
LD = struct with fields:
data_par1: [7680×14 double]
EEG_data = LD.data_par1(:,1);
Fs = 128;
Fn = Fs/2;
[th,df1] = bandpass(EEG_data, [4 7], Fs, 'ImpulseResponse','iir');
[fTheta,f2] = bandpass(EEG_data, [4 7], Fs, 'ImpulseResponse','iir');
[fAlpha,f3] = bandpass(EEG_data, [7 12], Fs, 'ImpulseResponse','iir');
[fBeta,f4] = bandpass(EEG_data, [12 25], Fs, 'ImpulseResponse','iir');
% calculate PSD
Theta = cpsd(fTheta,fTheta,5*Fs,2.5*Fs,5*Fs,Fs);
Alpha = cpsd(fAlpha,fAlpha,5*Fs,2.5*Fs,5*Fs,Fs);
Beta = cpsd(fBeta,fBeta,5*Fs,2.5*Fs,5*Fs,Fs);
% CPSD Plots
figure
subplot(3,1,1)
cpsd(fTheta,fTheta,5*Fs,2.5*Fs,5*Fs,Fs);
title('Theta')
subplot(3,1,2)
cpsd(fAlpha,fAlpha,5*Fs,2.5*Fs,5*Fs,Fs);
title('Alpha')
subplot(3,1,3)
cpsd(fBeta,fBeta,5*Fs,2.5*Fs,5*Fs,Fs);
title('Beta')
% calculate the average PSD of each signal
psdth = mean(Theta)
psdth = 1.2441
psda = mean(Alpha)
psda = 0.3879
psdb = mean(Beta)
psdb = 0.2075
Also, you are not actually doing any cross-signal calculations. Why not just use pwelch or one of the other options?
.

类别

在 Help Center 和 File Exchange 中查找有关 Pulsed Waveforms 的更多信息

标签

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by