www.gusucode.com > signal 工具箱matlab源码程序 > signal/medfreq.m

    function [f,pwr] = medfreq(varargin)
%MEDFREQ Median Frequency
%   FREQ = MEDFREQ(X) computes the median normalized angular frequency,
%   FREQ, of the power spectrum of the time-domain signal in vector X.
%   FREQ has units of radians/seconds.  If X is a matrix, MEDFREQ computes
%   the median frequency of each column in X independently.  MEDFREQ uses a
%   rectangular window when computing the spectrum.
%
%   The median frequency is defined as the frequency at which the power
%   spectrum is divided into two equal areas via rectangular integral
%   approximation.
%
%   FREQ = MEDFREQ(X, Fs) computes the median frequency, FREQ, of the power
%   spectrum of the time-domain signal in vector X with sample rate, Fs.
%   FREQ and Fs have units of hertz.
%
%   FREQ = MEDFREQ(Pxx, F) computes the median frequency of the power 
%   spectral density estimate, Pxx.  F is a vector containing the
%   frequencies that correspond to the estimates given in Pxx and must
%   contain at least two elements.
%
%   FREQ = MEDFREQ(Sxx, F, RBW) computes the median frequency of the power
%   spectrum estimate, Sxx, with resolution bandwidth RBW.
%
%   FREQ = MEDFREQ(..., FREQRANGE) specifies FREQRANGE as a two-element
%   vector of real values, specifying the two frequencies between which you
%   want to compute the median frequency.  The default value for FREQRANGE
%   is the entire bandwidth of the input signal.
%
%   [FREQ,PWR] = MEDFREQ(...) also returns the bandpower, POWER, of the
%   spectrum.  If FREQRANGE is specified, then POWER will contain the
%   bandpower within the frequency range.
%
%   MEDFREQ(...) with no output arguments will plot the PSD (or power
%   spectrum) and annotate the median frequency.
%
%   % Example 1:
%   %   Compute the median frequency of a chirp signal
%
%   nSamp = 1024;
%   Fs = 1024e3;
%   t = (0:nSamp-1)'/Fs;
%   x = chirp(t,50e3,nSamp/Fs,100e3);
%
%   medfreq(x,Fs)
%
%   % Example 2:
%   %   Compute the median frequency of a sinusoid from a PSD estimate
%
%   nSamp = 1024;
%   Fs = 1024e3;
%   t = (0:nSamp-1)'/Fs;
%   x = sin(2*pi*t*100.123e3);
%
%   [Pxx, F] = periodogram(x,kaiser(nSamp,38),[],Fs);
%   medfreq(Pxx,F)
%
%   See also MEANFREQ BANDPOWER FINDPEAKS PERIODOGRAM PWELCH PLOMB

%   Copyright 2014 The MathWorks, Inc.
narginchk(1,4);

% use a rectangular window for time-domain input
kaiserBeta = 0;

% fetch the PSD from the input
[Pxx, F, Frange, rbw, extraArgs, status] = psdparserange('medfreq', kaiserBeta, varargin{:});

% use full range if unspecified
if isempty(Frange)
  Frange = [F(1) F(end)];
end

% ensure no additional arguments are specified
if ~isempty(extraArgs)
  error(message('signal:medfreq:ExtraArgs'));
end

% compute the median frequency and power within the specified range
[f,pwr] = computeMedFreq(Pxx, F, Frange);

% plot if no output arguments specified
if nargout==0
  plotMedFreq(Pxx, F, Frange, rbw, f, status);
end

%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
function [f,pwr] = computeMedFreq(Pxx, F, freqrange)

% Compute the power from the PSD
width = specfreqwidth(F);
P = bsxfun(@times,width,Pxx);

% Cumulative rectangular integration
cumPwr = [zeros(1,size(P,2)); cumsum(P)];

% place borders halfway between each estimate.
cumF = [F(1); (F(1:end-1)+F(2:end))/2; F(end)];

% find the integrated power for the low and high frequency range
Plo = interpPower(cumPwr,cumF,freqrange(1));
Phi = interpPower(cumPwr,cumF,freqrange(2));

% return the power between the frequency range
pwr = Phi-Plo;

% return the frequency that divides the power equally
f = interpFreq(cumPwr,cumF,(Plo+Phi)/2);

%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
function p = interpPower(cumPwr,cumF,f)
idx = find(f<=cumF, 1,'first');
if ~isempty(idx)
  if idx==1
    p = signal.internal.linterp(cumPwr(1,:),cumPwr(2,:),cumF(1),cumF(2),f);
  else
    p = signal.internal.linterp(cumPwr(idx,:),cumPwr(idx-1,:), ...
                                cumF(idx),cumF(idx-1),f);
  end
else
  p = nan(1,size(cumPwr,2));
end

%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
function f = interpFreq(cumPwr, cumF, pwrThresh)
nChan = size(cumPwr,2);
f = zeros(1,nChan);

for iChan=1:nChan
  idx = find(pwrThresh(iChan)<=cumPwr(:,iChan),1,'first');
  if ~isempty(idx)
    if idx==1
       idx=2;
    end
    f(iChan) = signal.internal.linterp(cumF(idx-1), cumF(idx), ...
                 cumPwr(idx-1,iChan), cumPwr(idx,iChan), pwrThresh(iChan));
  else
    f(iChan) = NaN;
  end
end

%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
function plotMedFreq(Pxx, F, Frange, rbw, Fmed, status)

% plot spectrum
if strcmp(status.inputType,'power');
  % power spectrum when specified
  [hLine, xscale] = psdplot(Pxx, F, rbw, 'power', status);
else
  % otherwise, default to PSD
  [hLine, xscale] = psdplot(Pxx, F, rbw, 'psd', status);
end

% show the active frequency range of the measurement
hAxes = ancestor(hLine(1),'axes');
xLim = [F(1) F(end)];
yLim = get(hAxes,'YLim');
psdmaskactiverange(hAxes, xscale, xLim, yLim, Frange);

% plot vertical bar for each estimate
for i=1:numel(Fmed)
  line(xscale*[Fmed(i) Fmed(i)], yLim, ...
       'Parent',hAxes, ...
       'LineStyle','-.', ...
       'Color',get(hLine(i),'Color'));
end

% title the plot
titleStr = getString(message('signal:medfreq:MedianFreqEstimate'));
if isscalar(Fmed)
  [Fm, ~, units] = engunits(Fmed(1), 'unicode');
  if status.normF
    titleStr = sprintf('%s: %.3f \\times \\pi %srad/sample',titleStr,Fm/pi,units);
  else
    titleStr = sprintf('%s: %.3f %sHz',titleStr,Fm,units);
  end
end

title(titleStr);