Find the max psd

26 views (last 30 days)
vo
vo on 9 Dec 2022
Answered: William Rose on 10 Dec 2022
I have file data eeg after band pass contain one row.What code i can use to find max psd of this data.
  3 Comments
vo
vo on 9 Dec 2022
Edited: vo on 9 Dec 2022
Thank you, I run and it work but i don't know what exactly 5*fs 2.5*fs mean.
[pyy,f] = cpsd(y,y,5*fs,2.5*fs,5*fs,fs);
Can i change it?
William Rose
William Rose on 10 Dec 2022
@vo, you are welcome.
The line
[pyy,f] = cpsd(y,y,5*fs,2.5*fs,5*fs,fs);
estimates the power spectral density by Welch's method. Please read the help for cpsd() for more information. Yes, you can change those parameters.
The first "5*fs" argument tells cpsd to use a window of length 5*fs elements. BY default, cpsd uses a Hamming window. SInce fs is the sampling rate, a window of width 5*fs points corresponds to a 5 second long window. This means the frequency resolution of the PSD will be 1/(5 seconds) = 0.2 Hz. In other words, the spacing of the points in the PSD will be 0.2 Hz. If you increase the window width, the frequency spacing will be narrower. For example, if you select 10*fs (i.e. 10 second window width), the frequency resolution of the PSD will be 0.1 Hz.
The second argument, 2.5*fs, specifies the number of points to offset successive windows. This corresponds to a 2.5 second long offset between windows. This means each Hamming window is half-overlapped with the window before and the window after. Since the total signal duraiotn is 60 seconds, this means there wil be 12 5-second long windows that do not overlap, and another 11 in between those that half overlap, for a total of 23 windows. A PSD estimate will be generated form each window, and the 23 estimates will be combined to produce an estimate that has greater statistical accuracy than what you would get if you just had a single 60-second-long window.
You can change the 5*fs or the 2.5*fs and observe how the spectrum estimate changes.

Sign in to comment.

Accepted Answer

William Rose
William Rose on 10 Dec 2022
[This answer is a copy of my comment above. I intended my comment to be answer, but I posted it in the wrong place.]
I assume from your comment that you have the EEG for one channel, and it is a row vector. Let us assume that it has been sampled at 240 Hz for one minute.
First you must compute the power spectrum. Then you find the maximum power in the spectrum, and the associated frequency.
Example (with simulated data):
%Create a simulated EEG signal
fs=240; %sampling frequency (Hz)
td=60; %duration (s)
N=td*fs; %number of samples
y=randn(1,N); %Gaussian white noise
fhi=30; flo=0.5; %frequencies for bandpass filter (Hz)
[b,a]=butter(2,[flo/(fs/2),fhi/(fs/2)],'bandpass'); %create bandpass filter
y=filtfilt(b,a,y); %bandpass filter to make it more like an EEG
%Compute the power spectral density of the EEG
[pyy,f] = cpsd(y,y,5*fs,2.5*fs,5*fs,fs);
%Find the max power and the associated frequency
[pmax,idx]=max(pyy); %max power
fmax=f(idx); %frequency at max power
%Plot results
figure;
plot(f,pyy,'-b',fmax,pmax,'*r');
grid on; xlabel('Frequency (Hz)'); ylabel('PSD')
legend('PSD','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.=0.0121 at 5.6 Hz.
Try it. Good luck.

More Answers (0)

Categories

Find more on EEG/MEG/ECoG in Help Center and File Exchange

Tags

Community Treasure Hunt

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

Start Hunting!