Median Frequency

I am trying to use Matlab to find the median frequency of an EMG signal. The code I am using is below, and s is my signal. X is the frequency value I am trying to find so Z=P. Is there any solve function I can use to quickly find this answer or is guess and check my best bet? Or does anyone know an easieer method to find the median freq.?

>> Fs = 10000;

>> nfft = 2*nextpow2 (length(s));

>> Pxx = abs(fft(s,nfft)).^2/length(s)/Fs;

>> Hpsd = dspdata.psd(Pxx(1:length(Pxx)/2),'Fs',Fs);

>> plot(Hpsd);

>> P = avgpower(Hpsd)/2

>> Z = avgpower(Hpsd,[0,X])

Thanks

Answer 1

[Honglei Chen]

If I understand correctly, you want to find the frequency below which you have 50 percent of the signal power. If that's the case, you can use spectrum object to calculate the periodogram and then use cumsum to get the cumulative power distribution. You then just need to find the frequency corresponding to the 50th percetile of the power distribution, using your parameters as an example,

h = spectrum.periodogram;
Hpsd = psd(h,s,'Fs',1000,'NFFT',2*nextpow2(length(s)));
Pdist = cumsum(Hpsd.Data);
Freq = Hpsd.Frequencies;
OverHalfIdx = find(Pdist>=Pdist(end)/2,1,'first');
UnderHalfIdx = find(Pdist<=Pdist(end)/2,1,'last');
MidFreq = (Freq(OverHalfIdx)+Freq(UnderHalfIdx))/2

HTH

[da2841]

The results I have after executing this is,

OverHalfIdx = 1

UnderHalfIdx = Empty matrix: 0-by-1

MidFreq = Empty matrix: 0-by-1

But you are right in that I am trying to find the frequency below which I have 50% of the signal power.

Answer 2

[Wayne King]

I suspect your problem is that your signal has a nonzero mean possibly as a result of some nonstationarity or simply due to a DC shift. This would cause the OverHalfIdx variable in Honglei's code to be 1 and UnderHalfIdx to be empty.

I suggest you detrend your signal first. If you signal exhibits simply a DC shift, this can be accomplished with:

s = detrend(s,0);

If the nonzero mean is due to nonstationarity, it is trickier to handle.

Honglei's suggestion is a great solution after that.

[da2841]

Thanks Wayne and Honglei!

Answer 3

[mohamed ashour taha]

how can i use it in for loop % Hpsd = psd(h,S2,'Fs',Fs,'NFFT',NFFT);