www.gusucode.com > signal 工具箱matlab源码程序 > signal/obw.m
function [bw, flo, fhi, pwr] = obw(varargin) %OBW Occupied bandwidth. % BW = OBW(X) computes the 99% occupied bandwidth, BW, of the input % signal vector, X. BW has units of radians/sample. If X is a matrix, % then OBW computes the bandwidth over each column in X independently. % % To compute the occupied bandwidth, OBW first internally computes a % power spectral density estimate, Pxx, using PERIODOGRAM and a Kaiser % window. Next, the PSD is integrated with a rectangular approximation. % The bandwidth is computed from the frequency intercepts where the % integrated power crosses 0.5% and 99.5% of the total power in the % spectrum. % % BW = OBW(X, Fs) returns the occupied bandwidth, BW, in units of % hertz. Specify the sample rate of the signal, Fs, as a positive real % scalar. % % BW = OBW(Pxx, F) computes the occupied bandwidth of the PSD estimate, % Pxx. F is a vector of frequencies that corresponds to the vector of % Pxx estimates. If Pxx is a matrix, then OBW computes the bandwidth % over each column in Pxx independently. % % BW = OBW(Sxx, F, RBW) computes occupied bandwidth of the power spectrum % estimate, Sxx. F is a vector of frequencies that corresponds to the % vector of Sxx estimates. If Sxx is a matrix, then OBW computes the % bandwidth over each column in Sxx independently. % % BW = OBW(...,FREQRANGE,P) specifies the percentage, P, of the total % signal power present in the occupied band. Specify P as a positive % scalar value less than 100. The bandwidth is computed from the % frequency intercepts where the integrated power crosses the (100-P)/2 % and (100+P)/2 percentages of the total power in the spectrum. % % [BW, Flo, Fhi] = OBW(...) also returns the left and right % frequency borders of the occupied bandwidth. % % [BW, Flo, Fhi, POWER] = OBW(...) also returns the power within the % occupied bandwidth, POWER. % % OBW(...) with no output arguments by default plots the PSD (or power % spectrum) in the current figure window and annotates the bandwidth. % % % Example: % % Compute the occupied bandwidth of a chirp signal % % nSamp = 1024; % Fs = 1024e3; % t = (0:nSamp-1)'/Fs; % x = chirp(t,50e3,nSamp/Fs,100e3); % % obw(x,Fs) % % See also MEDFREQ POWERBW BANDPOWER FINDPEAKS PERIODOGRAM PWELCH PLOMB. % Copyright 2014 The MathWorks, Inc. narginchk(1,5); % use a rectangular window for time-domain input kaiserBeta = 0; % fetch the PSD from the input [Pxx, F, Frange, rbw, extraArgs, status] = psdparserange('obw', kaiserBeta, varargin{:}); % use full range if unspecified if isempty(Frange) Frange = [F(1) F(end)]; end % check if a percentage is specified if isempty(extraArgs) P = 99; else P = extraArgs{1}; validateattributes(P,{'numeric'},{'real','positive','scalar','<',100}, ... 'obw','P'); if numel(extraArgs)>1 error(message('signal:obw:UnrecognizedAdditionalArguments')); end end % compute the median frequency and power within the specified range [bw,flo,fhi,pwr] = computeOBW(Pxx, F, Frange, P); % plot if no output arguments specified if nargout==0 plotOBW(Pxx, F, Frange, rbw, flo, fhi, P, status); end %~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ function [bw,flo,fhi,pwr] = computeOBW(Pxx, F, Frange, P) % Compute the power from the PSD width = specfreqwidth(F); Pwr = bsxfun(@times,width,Pxx); % Cumulative rectangular integration cumPwr = [zeros(1,size(Pwr,2)); cumsum(Pwr)]; % 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,Frange(1)); Phi = interpPower(cumPwr,cumF,Frange(2)); % return the power between the frequency range totPwr = Phi-Plo; % return the frequency which intercepts the borders of the occupied band flo = interpFreq(cumPwr,cumF,Plo+(100-P)/200*totPwr); fhi = interpFreq(cumPwr,cumF,Plo+(100+P)/200*totPwr); % return the occupied bandwidth and occupied bandpower bw = fhi - flo; pwr = P/100 * totPwr; %~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 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 plotOBW(Pxx, F, Frange, rbw, flo, fhi, P, 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 translucent patch for each estimate for i=1:numel(flo) xData = xscale*[flo(i) fhi(i) fhi(i) flo(i)]; yData = yLim([1 1 2 2]); color = get(hLine(i),'Color'); patch(xData, yData, color, ... 'Parent',hAxes, ... 'EdgeColor','none', ... 'FaceAlpha',0.125); end % once patches are done, plot the frequency borders on top for i=1:numel(flo) line(xscale*[flo(i) flo(i)], yLim, ... 'Parent',hAxes, ... 'Color',get(hLine(i),'Color')); end for i=1:numel(fhi) line(xscale*[fhi(i) fhi(i)], yLim, ... 'Parent',hAxes, ... 'Color',get(hLine(i),'Color')); end % title the plot titleStr = getString(message('signal:obw:PercentOccupiedBandwidth',sprintf('%g',P))); if numel(flo)==1 [bw, ~, units] = engunits(fhi-flo, 'unicode'); if status.normF titleStr = sprintf('%s: %.3f \\times \\pi %srad/sample',titleStr,bw/pi,units); else titleStr = sprintf('%s: %.3f %sHz',titleStr,bw,units); end end title(titleStr);