Median Frequency
11 views (last 30 days)
Show older comments
da2841
on 23 Aug 2011
Answered: mohamed ashour taha
on 21 Dec 2015
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
0 Comments
Accepted Answer
Honglei Chen
on 23 Aug 2011
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
More Answers (2)
Wayne King
on 23 Aug 2011
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.
mohamed ashour taha
on 21 Dec 2015
how can i use it in for loop % Hpsd = psd(h,S2,'Fs',Fs,'NFFT',NFFT);
0 Comments
See Also
Categories
Find more on Spectral Measurements in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!