How to read spectrogram outputs and pick out interesting time-window?
22 views (last 30 days)
I am trying to determine where the strongest EEG activity happens, and how long it lasts (or what is the shortest time of this high activity which is enough to be recognized while implementing in BCI system for human use). I need to pick out the strongest activity and its lasting time, to save in a vector which I use later in classification of different kind of signals. For this I want to use PSD.
While reading literature review I noticed that many researchers uses Matlab’s spectrogram function and it looks like good one because of possibility to handle input parameters interesting for analyzing my EEG data. https://se.mathworks.com/help/signal/ref/spectrogram.html#bultmx7-s
I have attached EEG raw data that is 20 seconds long. The rest of the parameters are:
xy is the EEGdata attached to this question.
NFFT= length((xy)); %10240
If I run the spectrogram function with no output parameter I get the PSD of the signal presented in the picture below the code
It is easy to just zoom in and inspect visually where the activity is strong, but how do I do that in Matlab programmatically? I understand that spectrogram function gives PSD as well as a result in output parameters. So I have tried with this combination of outputs:
[s, f, t, ps] = spectrogram(xy,window,noverlap,NFFT,Fs,'yaxis' );
But the problem is that I need help to understand what the result is? What are rows and columns? Why are number of rows 5121 and columns 77? How can I use that information to get out the best activity and how long that activity lasts? If I understand it correctly, the s is returning complex sfft of the signal, but why is there 77 columns? I guess it has to do with windowing and the result of segmenting the signal, but how can I be sure that it is like that?
If I use s_real =real(s); I will get real part of fft, still 77 columns and 5121 rows.
If I use s_mag= abs(s); I will get magnitude of fft, still 77 columns and 5121 rows.
If I use PSD= 1/(Fs*10240).*s_mag.^2 I will get PSD, and still 77 columns and 5121 rows.
From question https://se.mathworks.com/matlabcentral/answers/280260-why-are-the-results-of-my-spectrogram-sinusoidal
In the description of spectrogram function, it stands that s matrix contains time and frequency: “time increases across the columns of s and frequency increases down the rows”. But the time in my resulting columns is not increasing. So what does the rows and columns representing? How can I use that data to find strongest activity and time window?
The output matrix ps from spectrogram function does not have the same result values as PSD result from PSD= 1/(Fs*10240).*s_mag.^2. Why are they different?
One more thing, if using spectrogram function with output parameters, the plot is empty- why is that? How can I plot PSD of frequency and time on x-axis? How can I plot chosen EEG activity PSD and its time using Matlab? If I use the code plot(t(:,:),PSD); I get following figure. What is represented on y-axis? I understand that I put PSD in the plot function, but how can I read from the figure what frequency give strongest value?
Later, when implementing chosen time window, for example, if I chose 1 second time window, how do I write the code to check for 2 second long time window? What I want with this is to inspect if the 2 second time window gives more accurate result then 1 second time window.
Any help will be gratefull!
LO on 3 Jul 2020
Edited: LO on 4 Jul 2020
the spectrogram can be handled as a surf plot
if f is the freq you get out of the spectrogram fs is the freq range of you segment of interest, p is the power from the spectrogram, Ref is the reference freq chosen for your segment, ts is a logic array calculated from your time vector (example if your spectrogram is over 60 sec time, your ts vector could be
[spec,f,t,p] = spectrogram(DATA,blackmanharris(nfft),overlap,nfft,Fs,'yaxis','MinThreshold',-80); % Fs is sampling freq. threshold value can be changed
ts = t > t(start) & t < t(end); % start and end are the indexes of the time vector your want to set as limits
fs = f > fs-50 & f< fs +300; % let's assume you want to select only a narrow band in your freq domain, you can also select all frequencies of course
s1.EdgeColor = 'none';
aspec = abs(spec); % remove the imaginary part of the spectrogram data
imagesc( t, f, log(1+aspec) );