what is the median frequency matlab code?

hello everybody,, my question what is the median frequency(MF) mat-lab code ,if the median frequency defined as the frequency that divides the magnitude-spectrum in two parts of equal size(the area under the curve for lower frequencies than MF equals the area under the curve for the frequencies higher than MF. first of all i've got the plot of the magnitude spectrum then i used trapz(x,y) to get the whole area under the curve, now i need a code to know what is the frequency that divides the total area into 2 equal parts ?!

[haya alawneh]

thanks for answering my question :),, i will try this code but actually my mission is to split EMG signal into 1 second for each window then to calculate the median frequency, relative power ,coherence , RMS and zero crossing rate where the signal change it's singe for each window, i used for loop to get 1 second frames , the second step is to implement the above code inside the for loop ? >> frame_duration=1; >> frame_len=frame_duration*fs; >> n=length(val); >> num_framed=floor(n/frame_len); >> fny=fs/2; >> for k=1:num_frames frameval((k-1)*frame_len + 1 :frame_len*k); frame val((k-1)*frame_len + 1 : frame_len*k); x_mags=abs(fft(val)); bin_vals=[0 : n-1]; fax_hz=bin_vals*fs/2; plot(fax_hz , x_mags) does it works ?!

[Star Strider]

Your unformatted code makes absolutely no sense to me.

I gave you code that will calculate the median frequency, as you requested.

[haya alawneh]

thankyou very much i tried the code :) ,, can i ask you another question s ?

[Star Strider]

My pleasure.

If my Answer solved your problem, please Accept it.

What is your question?

[haya alawneh]

the previous plot combine the frequency and the FFT amplitudes, now i need to square the FFT to get frquency-power plot then i need to calculate the relative power which is defined as the summation of the power in the 100-500 frequency range divided by the total power in the whole frequency domain? the second question is related to the zero crossing and how to check if the signal pass from positive to negative or from negative to positive ? can i use for loop to test the signals and save the index of zero crossing? Thankyou again

[Star Strider]

For relative power, use the find function to find the indices corresponding to your frequencies of interest, then use trapz to integrate the 100-500 range and the total power.

To calculate the power spectral density, see the documentation on Power Spectral Density Estimates Using FFT.

I do not know what you mean by ’zero crossing’. There should be no zero-crossing in the power spectral density.

[haya alawneh]

the zero crossing will be calculated from the amplitudes plot not the power plote ,,, can i use symsum to calculate the power? or should i use trapz function ?

[Star Strider]

I am still not certain what you are doing the zero-crossings on. You can use this little function to get their approximate indices:

zci = @(v) find(v(:).*circshift(v(:), [-1 0]) <= 0);                    % Returns Zero-Crossing Indices Of Argument Vector

Use trapz. (The cumtrapz function is also an option.)

[haya alawneh]

thankyou,, i will try to do this and sorry for many questions

[Star Strider]

My pleasure.

I don’t mind answering questions in my areas of expertise!

[haya alawneh]

according to the first code what is the difference between the time vector and the signal vector ?!

[Star Strider]

In the code in my Answer, the time vector (actually, deriving the sampling interval from it) is used to calculate the sampling frequency, and from the sampling frequency, the frequency vector.

The sampling interval, ‘Ts’ is defined in my code as:

Ts = mean(diff(t));

and the sampling frequency as:

Fs = 1/Ts;

I combined them in one statement in the code in my Answer.

Answer 1

[Star Strider]

You need to use the cumtrapz function:

t = ...;                                                % Time Vector
s = ...;                                                % Signal Vector

Fs = 1/mean(diff(t));                                   % Sampling Frequency
Fn = Fs/2;                                              % Nyquist Frequency
L = length(t);

FTs = fft(s)/L;
Fv = linspace(0, 1, fix(L/2)+1)*Fn;                     % Frequency Vector
Iv = 1:length(Fv);                                      % Index Vector

CumAmp = cumtrapz(Fv, abs(FTs(Iv)));                    % Integrate FFT Amplitude
MedFreq = interp1(CumAmp, Fv, CumAmp(end)/2);           % Use ‘interp1’ To Find ‘MF’

figure(1)
plot(Fv, abs(FTs(Iv))*2, '-b')                          % Plot FFT
hold on
plot(Fv, CumAmp, '-g')                                  % Plot Cumulative Amplitude Integral
plot([MedFreq MedFreq], ylim, '-r', 'LineWidth',1)      % Plot Median Frequency
hold off
grid

I tested this with an archive signal, so the code works. It should work with your signal without modification. You can delete the plotting of the cumulative amplitude integral (green line), since it is there only to illustrate that the median frequency is calculated correctly.

[Gonzalo Guerrero]

Hi,

This is a very good code, thank you for bring it up. However, when I use it, and compare with other software, my results are a bit off, most likely I am doing something very wrong. I will explain you the issue, maybe you can help and guide me.

I have a sinusoidal curve, that last for a few seconds and reach different values. When it comes to the analysis, I select my points of interes. So I ask the script to go point by point with a window of 500ms and a overlap of 499 ms.

Then, T for me is = to the 500ms (i.e: t=mypoints(nn:nn+962) )from the point that it is analysing. and S= the values of that time period too.

Could you help, please? :D

Thank you