Fix the code calculate PSD

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

[Tran]

data_par1(:,1) mean that i choose the 1 channels(AF3) to do

Answer 1

[Star Strider]

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?

.