Hi. I need some help about Spectrograms

I have EMG data, which I split into active contraction zones which I did and computed their spectrograms. How can I extract dominant and median frequency from the signal spectrum?

Apprectiate the help, thanks!

Answer 1

[Star Strider]

I would use the pspectrum function for this with the 'spectrogram' option, since the units and results are easier to interpret. It has several options for outputs, so choose the one that works best for what you want to do. I am not certain what frequencies you are referring to, or exactly what you want to do, however extracting the frequencies as a function of time is straightforward.

Example —

fs = 3000;
t = 0:1/fs:1-1/fs;
x2 = exp(2j*pi*100*cos(2*pi*2*t)) + randn(size(t))/100;

[p,f,t] = pspectrum(x2,fs,'spectrogram');

figure
surf(t,f,p,'EdgeColor','none')
colormap(turbo)
view(0,90)
xlabel('Time')
ylabel('Frequency')
zlabel('Power')

tr = t>0.009 & t<0.101;                         % Time Range
idx = find(tr);                                 % Indices
fv = f(idx)                                     % Frequency

fv = 5×1

1.0e+03 * -1.5000 -1.4971 -1.4941 -1.4912 -1.4883

pwr = p(:,idx);                                 % Power
maxpwr = max(pwr)                               % Maximum Power

maxpwr = 1×5

0.2213 0.2436 0.2854 0.3530 0.4771

.

[Kaosar Soomro]

Hey, thanks for your reply. I'm quite new to matlab so plesae bear with me!

So far this is what I have done:

1. I first split data into zones:

%for emg1_flexionextension

zone1_emg1_flexionextension = emg1_flexionextension (1420:2400);

zone2_emg1_flexionextension = emg1_flexionextension (2400:3500);

zone3_emg1_flexionextension = emg1_flexionextension (3800:4400);

zone4_emg1_flexionextension = emg1_flexionextension (5000:5500);

zone5_emg1_flexionextension = emg1_flexionextension (6100:6800);

zone6_emg1_flexionextension = emg1_flexionextension (6850:7350);

2. Then I computed spectrogram for eah of these zones: (frequency-time spectrogram)

% for emg1_flexionextension zones

fs = 1000;

dt = 1/fs;

t = [0:5000]*dt;

F25=figure;

spectrogram(zone1_emg1_flexionextension,rectwin(800),250,256,fs,'yaxis');

F26=figure;

spectrogram(zone2_emg1_flexionextension,rectwin(800),250,256,fs,'yaxis');

F27=figure;

spectrogram(zone3_emg1_flexionextension,rectwin(512),250,256,fs,'yaxis');

F28=figure;

spectrogram(zone4_emg1_flexionextension,rectwin(450),250,256,fs,'yaxis');

F29=figure;

spectrogram(zone5_emg1_flexionextension,rectwin(512),250,256,fs,'yaxis');

F30=figure;

spectrogram(zone6_emg1_flexionextension,rectwin(450),250,256,fs,'yaxis');

From these zones, how would I go about extracting dominant and median frequency?

Thanks for your help!

[Star Strider]

If you have the time values already determined, it may not be necessary to calculate the spectrogram. Just use the fft function and use the abs values of it.

I assume the dominant frquency would be the maximum frequency, and you can determine that and its index with the max function, and then compute the median frequency with the medfreq function.

[Kaosar Soomro]

OH right, for the data, its a continous EMG recording so every second its takes a reading. So for each contraction of muscle, it shows in the signal and for each of these contractions I was required to compute a spectrogram for each contraction zone and then extract the dominant and median frequecy from these zones.

Oh so it would be max (spectrogram)?

[Star Strider]

A sampling frequency of 1 Hz is likely too slow to be of any value with respect to analysing any physiological signal.

As I mentioned, if you are analysing the data at each time window, you are already calculating the Fourier transform of it so I doubt that the spectrogram is necessary. Just calculate the fft and medfreq of each time segment.

[Kaosar Soomro]

okay, this is great although trying to extract the median frequency from one zone in the EMG data - would I still go about it this way, or can i use the medfreq function somehow?

Would you mind breaking down what id need to input in each line as im unsure what most of these functions mean & what they do.

Thank you for your patience!

[Star Strider]

My pleasure!

I would not use spectrogram or anything like it unless it is absolutely necessary. However I doubt that it is necessary here.

If you have specific signal segments defined by time, calculate the fft of each segment and save the max frequency and medfreq on each segment. Those should be scalar values, so just save them in an array that you can reference with a spedific time index. The maximum frequency would be the frequency associated with the index of the maximum amplitude, so:

Fs = ...; % Sampling Frequency

Fn = Fs/2; % Nyquist Frequency

L = length(signal);

NFFT = 2^nextpow2(L);

FTsignal = fft(signal-mean(signal),NFFT)/L;

Fv = linspace(0, 1, NFFT/2+1)*Fn; % Frequency Vector

Iv = 1:numel(Fv); % Index Vector

figure

plot(Fv, mag2db(abs(FTsignal(Iv))*2))

grid

[amp,idx] = max(abs(FTsignal(Iv)))

domfrq = Fv(idx) % Dominant Frequency

medfrq = medfreq(signal,Fs); % Median Frequency

That should work for each chosen time segment, so it may be necessary to put it in a loop, especially since you want both the dominant and median frequencies.

.

[Kaosar Soomro]

Ah okay, it just says to "Compute the spectrogram of each active zone, and then extract the dominant frequency and the median frequency of the signal spectrum".

Yes, I have specific signal segments defined by time (pasted earlier).

Just a question, in these two lines,

Fv = linspace(0, 1, NFFT/2+1)*Fn;

Iv = 1:numel(Fv);

what do the numbers in (0, 1, NFFT/2+1) mean? and what does 1:numel mean? - just wondering just to make sure if its ok to select any values or if there has to be a specific number.

Thank you!

[Star Strider]

I am not certain what you are doing. I got the impression this was research.

‘what do the numbers in (0, 1, NFFT/2+1) mean?’

That line computes the frequency vector for the Fourier transform plot. It creates a vector going from 0 to 1 with‘NFFT/2+1’ elements, and then multiplies it by the Nyquist frequency.

The second ‘Iv’ vector simply indexes into the Fourier transfomr result.

.

[Kaosar Soomro]

Oh I see, is there a way to carry this out for each chosen time segment without putting it in a loop?

Alternatively, could I use [s,f,t] = spectrogram(x, Fs) to find dominant and median frequency?

Thank you

[Star Strider]

Considering that you already created the time ranges, the loop is likely the most efficient option, especially since you want the median frequency as well as the dominant frequency.

If you need the spectrogram, simply calculate it separately.

[Kaosar Soomro]

Ah okay, would this be by using the functions above?

Would you mind breaking down how to create a loop in order to use the function above?

Thank you

[Star Strider]

Your other post about segmenting the signal by time is part of it. Just calculate the fft of whatever collumns you want, and get the frequency of the maximum value for the dominant frequency and calculate the median frequency with medfreq separately for each segment. It is likely possible to do that with the entire buffer matrix at once with each function, and then either plot all those results or only the ones you are interested in.

Example —

EMG = randn(1,5E+5);                                                % Create Data
Fs = 1E+3;                                                          % Sampling Frequency
rowlen = 20*Fs;                                                     % Seconds * Samples/Second = Samples
EMGbuf = buffer(EMG,rowlen);
Sec40 = EMGbuf(:,fix(40/20));
Sec80 = EMGbuf(:,fix(80/20));
Sec120 = EMGbuf(:,fix(120/20));
t = linspace(0, rowlen-1, rowlen)/Fs;                               % Common Time Vector

Fn = Fs/2;                                                          % Nyquist Frequency
L = size(EMGbuf,1);
NFFT = 2^nextpow2(L);
FTsignal = fft(EMGbuf-mean(EMGbuf),NFFT)/L;
Fv = linspace(0, 1, NFFT/2+1)*Fn;                                   % Frequency Vector
Iv = 1:numel(Fv);                                                   % Index Vector

% figure
% plot(Fv, mag2db(abs(FTsignal(Iv,:))*2))
% grid
% ylim([-100 max(ylim)])

[amp,idx] = max(abs(FTsignal(Iv,:)));
domfrq = Fv(idx);                                                   % Dominant Frequency
medfrq = medfreq(EMGbuf,Fs);                                        % Median Frequency
buftime = linspace(1, size(EMGbuf,2), size(EMGbuf,2))*rowlen/Fs;

figure
yyaxis left
stairs(buftime, domfrq, 'DisplayName','Dominant Frequency')
ylabel('Frequency (Hz)')
yyaxis right
stairs(buftime, medfrq, 'DisplayName','Median Frequency')
grid
xlabel('Buffer Column (s)')
ylabel('Frequency (Hz)')
legend('Location','best')

This is how I would do it.

Use your own data vector for my ‘EMG’ vector.

.

[Star Strider]

I already did, as I illustrated. No spectrogram function that I can find (spectrogram, pspectrum) allows the median frequency as an option. If you want both, my code is the only possible approach.

[Star Strider]

Kaosar Soomro’s Answer moved here —

Ah okay, could I use the below functions on each segment (active zone where muscle contraction occurs)? - how would I know what to put for the timeColumn?

[s,f,t] = spectrogram(x, Fs);

timeColumn = 5;

intpwr = cumtrapz(f, abs(s(:,timeColumn)));

medianpwrfreq = interp1(intpwr, f, max(intpwr)/2);

Appreciate the help!

[Star Strider]

I am lost.

This is not compatible with the buffer matrix you want to use.

I have no idea how to make spectrogram or pspectrum compatible with the buffer matrix, other than to compute the spectrogram on each column of the buffer matrix. I demonstrated that elsewhere.

EDIT — (3 Jum 2022 at 16:12)

I can get specified time bins with pspectrum however not with spectrogram —

EMG = randn(1,2.5E+5);                                              % Create Data
Fs = 1E+3;                                                          % Sampling Frequency
L = numel(EMG);
tv = linspace(0, L-1, L)/Fs;

[sp,sf,st] = pspectrum(EMG, tv, 'spectrogram', 'TimeResolution',20, 'OverlapPercent',0);

figure
surfc(st,sf,sp, 'EdgeColor','none')
grid on
xlabel('Time (s)')
ylabel('Frequency (Hz)')
colormap(turbo)
view(0,90)
set(gca,'XTick',st+10, 'XTickLabel',st)                             % Ticks At Time Bin Midpoints, So 10, 30, ...

Neither one has a median frequency option.

.

Answer 2

[Kaosar Soomro]

This is great, However for each segment (these are not the same as in my other post) which I have defined as certain regions of the signal, I have generated a spectrogram (time-frequency), I'm interested in extracting the dominat and median frequency from each of these spectrograms.

Could you break down how would I go about doing this?

Thanks

Answer 3

[Kaosar Soomro]

Ah okay, could I use the below functions on each segment (active zone where muscle contraction occurs)? - how would I know what to put for the timeColumn?

[s,f,t] = spectrogram(x, Fs);

timeColumn = 5;

intpwr = cumtrapz(f, abs(s(:,timeColumn)));

medianpwrfreq = interp1(intpwr, f, max(intpwr)/2);

Appreciate the help!