www.gusucode.com > signal 工具箱matlab源码程序 > signal/+signal/+internal/getToneFromPSD.m
function [power, freq, idxTone, idxLeft, idxRight] = getToneFromPSD(Pxx, F, rbw, toneFreq) %GETTONEFROMPSD Retrieve the power and frequency of a windowed sinusoid % % This function is for internal use only and may be removed in a future % release of MATLAB % Copyright 2013 The MathWorks, Inc. if nargin<4 [~, idxTone] = max(Pxx); elseif F(1) <= toneFreq && toneFreq <= F(end) % find closest bin to specified freq [~, idxTone] = min(abs(F-toneFreq)); % look for local peak in vicinity of tone iLeftBin = max(1,idxTone-1); iRightBin = min(idxTone+1,numel(Pxx)); [~, idxMax] = max(Pxx(iLeftBin:iRightBin)); idxTone = iLeftBin+idxMax-1; else power = NaN; freq = NaN; idxTone = []; idxLeft = []; idxRight = []; return end % sidelobes treated as noise idxLeft = idxTone - 1; idxRight = idxTone + 1; % roll down slope to left while idxLeft > 0 && Pxx(idxLeft) <= Pxx(idxLeft+1) idxLeft = idxLeft - 1; end % roll down slope to right while idxRight <= numel(Pxx) && Pxx(idxRight-1) >= Pxx(idxRight) idxRight = idxRight + 1; end % provide indices to the tone border (inclusive) idxLeft = idxLeft+1; idxRight = idxRight-1; % compute the central moment in the neighborhood of the peak Ffund = F(idxLeft:idxRight); Sfund = Pxx(idxLeft:idxRight); freq = dot(Ffund, Sfund) ./ sum(Sfund); % report back the integrated power in this band if idxLeft<idxRight % more than one bin power = bandpower(Pxx(idxLeft:idxRight),F(idxLeft:idxRight),'psd'); elseif 1 < idxRight && idxRight < numel(Pxx) % otherwise just use the current bin power = Pxx(idxRight) * (F(idxRight+1) - F(idxRight-1))/2; else % otherwise just use the average bin width power = Pxx(idxRight) * mean(diff(F)); end % protect against nearby tone invading the window kernel if nargin>2 && power < rbw*Pxx(idxTone) power = rbw*Pxx(idxTone); freq = F(idxTone); end