@vo,
[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:
y1=cell2mat(struct2cell(load('EEGdata')));
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)
plot(f,pyy,'-b',fmax,pmax,'*r',fmaxA,pmaxA,'*g');
grid on; xlabel('Frequency (Hz)'); ylabel('PSD')
legend('PSD','Max','\alpha max'); title('EEG Power Spectral Density')
subplot(212)
plot(f,pyy,'-b',fmax,pmax,'*r',fmaxA,pmaxA,'*g');
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)
fprintf('Alpha band: max P.S.D.=%.4f at %.1f Hz.\n',pmaxA,fmaxA)
Try it.
