www.gusucode.com > signal 工具箱matlab源码程序 > signal/+fdesign/@abstractmeas/findfrequency.m

    function F = findfrequency(this, hfilter, A, direction, place)
%FINDFREQUENCY   Find the frequency point for the given amplitude
%   F = FINDFREQUENCY(H, HD, A, DIR, PLACE) returns the frequency point F
%   for the amplitude A in the filter object HD.  A is in linear units and
%   is relative to the nominal gain of the filter.
%
%   F = FINDFREQUENCY(H, 


%   Copyright 2005 The MathWorks, Inc.

narginchk(4,5);

N = 2^10;

ngain = nominalgain(hfilter);
if ~isempty(ngain)
    A = ngain*A;
end

if this.NormalizedFrequency, Fs = 2;
else,                        Fs = this.Fs; end

if strcmpi(direction, 'up'), up = true;
else,                        up = false; end

if strcmpi(place, 'first'), first = true;
else,                       first = false; end

% find first point that fails a3db
[h,w] = freqz(hfilter, N, Fs);
h = abs(h);

[w_lo, w_hi] = lclfindfrequency(w, h, A, up, first);

[F htest] = getF(hfilter, [w_lo w_hi], Fs);

count = 1;

while abs(A-htest) > 1e-5 && count < 4

    % Do a second refinement if we did not hit the point exactly.

    [h, w] = freqz(hfilter, linspace(w_lo, w_hi, N), Fs);
    h = abs(h);

    [w_lo, w_hi] = lclfindfrequency(w, h, A, up, first);

    [F htest] = getF(hfilter, [w_lo w_hi], Fs);
    
    count = count+1;
end

% -------------------------------------------------------------------------
function [w_lo, w_hi] = lclfindfrequency(w, h, A, up, first)

if first
    if up
        indx = find(h >= A, 1, 'first');
    else
        indx = find(h <= A, 1, 'first');
    end
    if isempty(indx), indx = length(w); end

    w_lo = w(max(1, indx-1));
    w_hi = w(indx);
else
    if up
        indx = find(h <= A, 1, 'last');
    else
        indx = find(h >= A, 1, 'last');
    end
    if isempty(indx), indx = 1; end
    w_lo = w(indx);
    w_hi = w(min(length(w), indx+1));
end


% -------------------------------------------------------------------------
function [F, htest] = getF(hfilter, w, Fs)

F = mean(w);

htest = freqz(hfilter, [0 F], Fs);
htest = abs(htest(2));

% [EOF]