PSD with Alpha(8-13HZ)

49 views (last 30 days)
vo on 12 Dec 2022
Edited: William Rose on 12 Dec 2022
I want to caculate the max PSD of Alpha wave of EEG signal. What problem with my code(PSD too high) ?
In this case, I have a raw data eeg in 1 channel (I show file data below) record at fs =128 and this is my code.
I think my mistake in pwelch function, Anyone can help me?
% i want band pass 2 time ( 4-30) --> (8-12Hz)
Band4_30 = bandpass(Data(1,:), [4 30], 128, 'ImpulseResponse','iir');
Alpha = bandpass(Band4_30(1,:), [8 12], 128);
%Compute the power spectral density of the EEG
%[pb,f] = pwelch(Alpha(1,:),[],[],[],128);
[pb,f] = pwelch(Alpha(1,:),hamming(256),[],256,128);
%Find the max power
%Display results
fprintf('Maximum PSD of Alpha=%.4f',pmax)

Answers (1)

William Rose
William Rose on 12 Dec 2022
Edited: William Rose on 12 Dec 2022
[edit: replace 8 with flo and 12 with fhi in the code to compute vector fA.]
The data file you gave is a structure. The first line loads it as a regular numeric array, because I find it easier to work with a regular array.
To find the max power point in the alpha wave frequency band (8-12 Hz), compute the PSD using the code I gave in repsonse to your other question, then find the max in the 8-12 Hz band. The code fragment below assumes that pyy is the vector containing the P.S.D. estimate, and it assumes f is the vector of frequencies associated with pyy.
The code below generates two figures. Figure 1 shows the EEG versus time, before and after highpass filtering. Figure 2 shows the PSD, and a zoomed-in view of the PSD. The zoomed-in plot shows that there is very little power in the alpha band. compared to the power at lower frequencies.
% flo=8; fhi=12; %frequency limits (Hz)
% [pmaxA,idxA]=max(pyy(f>=flo & f<=fhi)); %max power in band
% fA=f(f>=8 & f<=12); %frequency vector, in band
% fmaxA=fA(idxA); %frequency of max power in band
The plot of the full length recording (below) shows a spike at about t=15 s, and it shows a non-zero and drifting baseline. The recording is dominated by a strong non-physiological oscillation at 28 Hz. The sampling rate is half of the minimum rate recommended by the American Clinical Neurophysiology Society (Guideline 4, 2016).
The recording is 199 seconds long. I will apply Welch's method with half-overlapped windows of length 10 seconds. This means the final 4 seconds of the recording will not be used in the power spectrum computation, because the last window that fits will be from t=185 to t=195 seconds.
%Read EEG from file
%File is a struct, so convert it to a regular numeric array:
fs=128; %sampling rate (Hz)
t=(0:length(y1)-1)/fs; %time vector (s)
%Apply a highpass filter to the data to remove DC offset and baseline drift.
fco=0.5; %highpass filter cutoff frequency (Hz)
[b,a]=butter(2,fco/(fs/2),'high'); %create highpass filter
y=filtfilt(b,a,y1); %apply filter to raw signal
%Plot the signal and a zoomed-in segment and signal after filtering.
figure; subplot(311); plot(t,y1,'-b');
xlabel('Time (s)'); title('Raw Signal')
subplot(312); plot(t,y1,'-b');
xlim([100 102]); xlabel('Time (s)'); title('Raw Signal')
subplot(313); plot(t,y,'-b');
xlabel('Time (s)'); title('Highpass Filtered Signal')
%Compute the power spectral density of the EEG
[pyy,f] = cpsd(y,y,10*fs,5*fs,10*fs,fs);
%Find the max power and the associated frequency
[pmax,idx]=max(pyy); %max power
fmax=f(idx); %frequency at max power
%Find max power and associated frequency in alpha band (8-12 Hz)
flo=8; fhi=12; %frequency limits (Hz)
[pmaxA,idxA]=max(pyy(f>=flo & f<=fhi)); %max power in band
fA=f(f>=flo & f<=fhi); %frequency vector, in band
fmaxA=fA(idxA); %frequency of max power in band
%Plot results
figure; subplot(211)
grid on; xlabel('Frequency (Hz)'); ylabel('PSD')
legend('PSD','Max','\alpha max'); title('EEG Power Spectral Density')
grid on; xlabel('Frequency (Hz)'); ylabel('PSD'); xlim([0,20]);
legend('PSD','Max','\alpha max'); title('EEG Power Spectral Density')
%Display results on console
fprintf('Maximum EEG P.S.D.=%.4f at %.1f Hz.\n',pmax,fmax)
Maximum EEG P.S.D.=263.4882 at 27.9 Hz.
fprintf('Alpha band: max P.S.D.=%.4f at %.1f Hz.\n',pmaxA,fmaxA)
Alpha band: max P.S.D.=0.2901 at 8.4 Hz.
Try it.


Find more on EEG/MEG/ECoG 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!