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

    function pwr = bandpower(varargin)
%BANDPOWER Band power
%   P = BANDPOWER(X) computes the average power in the input signal vector,
%   X.  If X is a matrix, then BANDPOWER computes the average power in each
%   column independently.
%
%   P = BANDPOWER(X, Fs, FREQRANGE) specifies FREQRANGE as a two-element
%   vector of real values, specifying the two frequencies between which you
%   want to measure the power.  Fs is the sampling rate of the input
%   signal.  The power is estimated by applying a Hamming window and using
%   a periodogram of the same length as the input vector.  If the input
%   vector X contains N samples, then FREQRANGE must be contained within
%   the interval:
%      [              0, Fs/2          ] if X is real and N is even
%      [              0, Fs*(N-1)/(2*N)] if X is real and N is odd
%      [-Fs*(N-2)/(2*N), Fs/2          ] if X is complex and N is even
%      [-Fs*(N-1)/(2*N), Fs*(N-1)/(2*N)] if X is complex and N is odd
%
%   P = BANDPOWER(Pxx, F, 'psd') computes the average power via a rectangle
%   approximation of the integral of the Power Spectral Density (PSD)
%   estimate, given in the vector Pxx over the frequencies specified in
%   vector F.
%
%   P = BANDPOWER(Pxx, F, FREQRANGE, 'psd') specifies FREQRANGE as a
%   two-element vector of real values, specifying the two frequencies
%   between which you want to measure the power.  F is a vector containing
%   the frequencies that correspond to the estimates given in Pxx.
%
%   NOTE: If the frequency range values don't exactly match the frequency
%   values stored in F then the next closest value is used.
%
%   % Example
%   %   Examine the power present within 50 kHz of the carrier of a 1 kHz
%   %   sinusoid FM modulated at 250 kHz with a modulation index of 0.1.
%
%   Fs = 1.0e6; Fc = 250e3; Fd = 1e3; m = 0.1;
%   x = vco(m*sin(2*pi*Fd*(1:1000)/Fs),Fc,Fs);
%   periodogram(x,[],length(x),Fs);
%   mod_power_dB = 10*log10(bandpower(x, Fs, Fc + [-50e3 50e3]))
%
%   See also POWERBW OBW PERIODOGRAM PWELCH MEANFREQ MEDFREQ PLOMB.

%   Copyright 2012-2014 The MathWorks, Inc.

narginchk(1,4);

matches = find(strcmpi('psd',varargin));
varargin(matches) = [];

if any(matches)
    pwr = psdbandpower(varargin{:});
else
    pwr = timedomainbandpower(varargin{:});
end

function pwr = timedomainbandpower(x, fs, freqrange)

% perform column vector conversion before checking 2d matrix
if isvector(x)
  x = x(:);
end

validateattributes(x, {'numeric'},{'2d','finite'}, ...
    'bandpower', 'x', 1);
  
if nargin==1
    % full range specified
    pwr = rms(x).^2;
    return
elseif nargin == 2
    error(message('signal:bandpower:FreqRangeMissing'));
else
    validateattributes(fs, {'numeric'},{'scalar','finite','real','positive'}, ...
        'bandpower', 'Fs', 2);
    % Compute periodogram using a hamming window with the same length as
    % input
    
    % Cast to enforce Precision rules
    fs = double(fs);
    freqrange = signal.internal.sigcasttofloat(freqrange,'double',...
      'bandpower','FREQRANGE','allownumeric');
    
    n = size(x,1);
    
    if isreal(x)
        [Pxx, F] = periodogram(x, hamming(n), n, fs);
    else
        [Pxx, F] = periodogram(x, hamming(n), n, fs, 'centered');
    end
    % Return the bandpower
    pwr = psdbandpower(Pxx, F, freqrange);
end

function pwr = psdbandpower(Pxx, W, freqrange)

% perform column vector conversion before checking 2d matrix
if isvector(Pxx)
  Pxx = Pxx(:);
end

validateattributes(Pxx, {'numeric'},{'2d','finite','real'}, ...
    'bandpower', 'Pxx', 1);
  
if nargin < 2
    error(message('signal:bandpower:FreqVectorMissing'));
end

validateattributes(W,{'numeric'},{'vector','finite','real'}, ...
    'bandpower', 'F', 2);
if size(Pxx,1) ~= numel(W)
    error(message('signal:bandpower:FreqVectorMismatch'));
end

if any(diff(W)<0),
    error(message('signal:bandpower:FreqVectorMustBeStrictlyIncreasing'));
end

% Cast to enforce Precision rules
W = double(W);
if nargin < 3,        
    freqrange = [W(1) W(end)];
    freqrangespecified = false;
else
    validateattributes(freqrange,{'numeric'},{'vector','finite','real'}, ...
        'bandpower', 'FREQRANGE', 3);    
    if length(freqrange)~=2
        error(message('signal:bandpower:FreqRangeTwoElement'));
    elseif freqrange(1)<W(1) || freqrange(2)>W(end)
        error(message('signal:bandpower:FreqRangeOutOfBounds'));
    elseif freqrange(1) >= freqrange(2)
        error(message('signal:bandpower:FreqRangeMustBeStrictlyIncreasing'));
    end
    freqrangespecified = true;
end

% force column vector
W = W(:);

% Find indices of freq range requested.
idx1 = find(W<=freqrange(1), 1, 'last' );
idx2 = find(W>=freqrange(2), 1, 'first');

% Determine the width of the rectangle used to approximate the integral.
width = diff(W);
if freqrangespecified
    lastRectWidth = 0;  % Don't include last point of PSD data.
    width = [width; lastRectWidth];
else
    % There are two cases when spectrum is twosided, CenterDC or not.
    % In both cases, the frequency samples does not cover the entire
    % 2*pi (or Fs) region due to the periodicity.  Therefore, the
    % missing freq range has to be compensated in the integral.  The
    % missing freq range can be calculated as the difference between
    % 2*pi (or Fs) and the actual frequency vector span.  For example,
    % considering 1024 points over 2*pi, then frequency vector will be
    % [0 2*pi*(1-1/1024)], i.e., the missing freq range is 2*pi/1024.
    %
    % When CenterDC is true, if the number of points is even, the
    % Nyquist point (Fs/2) is exact, therefore, the missing range is at
    % the left side, i.e., the beginning of the vector.  If the number
    % of points is odd, then the missing freq range is at both ends.
    % However, due to the symmetry of the real signal spectrum, it can
    % still be considered as if it is missing at the beginning of the
    % vector.  Even when the spectrum is asymmetric, since the
    % approximation of the integral is close when NFFT is large,
    % putting it in the beginning of the vector is still ok.
    %
    % When CenterDC is false, the missing range is always at the end of
    % the frequency vector since the frequency always starts at 0.
    
    % assuming a relatively uniform interval
    missingWidth = (W(end) - W(1)) / (numel(W) - 1);
    
    % if CenterDC was not specified, the first frequency point will
    % be 0 (DC).
    centerDC = ~isequal(W(1),0);
    if centerDC
        width = [missingWidth; width];
    else
        width = [width; missingWidth];
    end
end

% Sum the average power over the range of interest.
pwr = width(idx1:idx2)'*Pxx(idx1:idx2,:);

% [EOF]