fdbutter.m

matlabDigitalSigalProcess内有文件若干

        fChangeMeas(i,type,Fs,Lresp,L1,L2,f1,f2,f3,f4,Rp,Rs)
        return
    end
    ax.xlimPassband = fdutil('xlimpassband',type,...
                         Fs,f1.value,f2.value,f3.value,f4.value);

    % fdutil('updateRanges',i,f1,f2,f3,f4)

%---------------------------------------------------------------------
% fdbutter('L2drag',type,ind)
%  drag callback of L2 - stopband line
%  Inputs:
%     type - band configuration 1==low, 2=high, 3=pass, 4=stop
%     ind - index of vertex being dragged
case 'L2drag'
    type = varargin{2};
    ind = varargin{3};
    xd = L2.xdata;
    minspacing = Fs/500;
    switch type
    case 1  % lowpass
        newf2 = inbounds(xd(ind),[f1.value+minspacing Fs/2-minspacing]);
        xd(1) = newf2;
        L2.xdata = xd;
        f2.value = newf2;
        i = 2;
    case 2  % highpass
        newf1 = inbounds(xd(ind),[minspacing f2.value-minspacing]);
        xd(2) = newf1;
        L2.xdata = xd;
        f1.value = newf1;
        i = 1;
    case 3  % bandpass
    % L2 xdata = [0 f1 NaN f4 Fs/2]
        if ind == 2   % dragging f1
            newf1 = inbounds(xd(ind),[minspacing f2.value-minspacing]);
            xd(2) = newf1;
            L2.xdata = xd;
            f1.value = newf1;
            i = 1;
        else % dragging f4
            newf4 = inbounds(xd(ind),[f3.value+minspacing Fs/2-minspacing]);
            xd(4) = newf4;
            L2.xdata = xd;
            f4.value = newf4;
            i = 4;
        end
    case 4   % bandstop
    % L2 xdata = [f2 f3]
        if ind == 1   % dragging f2
            newf2 = inbounds(xd(ind),[f1.value+minspacing f3.value-minspacing]);
            xd(1) = newf2;
            L2.xdata = xd;
            f2.value = newf2;
            i = 2;
        else % dragging f3
            newf3 = inbounds(xd(ind),[f2.value+minspacing f4.value-minspacing]);
            xd(2) = newf3;
            L2.xdata = xd;
            f3.value = newf3;
            i = 3;
        end
    end
    if ~minOrdCheckbox.value
        fChangeMeas(i,type,Fs,Lresp,L1,L2,f1,f2,f3,f4,Rp,Rs)
    end
    
%---------------------------------------------------------------------
% fdbutter('Rpdrag',type,ind)
%  drag callback of L1 - passband line
%  Inputs:
%     type - band configuration 1==low, 2=high, 3=pass, 4=stop
%     ind - index of segment being dragged
case 'Rpdrag'
    type = varargin{2};
    ind = varargin{3};
    yd = L1.ydata;
    
    below = yd(ind);
    if below >= 0 
        below = -.00001;
    elseif below < -10
        below = -10;
    end
    above = 0;
    
    newRp = above-below;

    if minOrdCheckbox.value
        if type ~= 4
        % 'ydata':[maxpass maxpass NaN minpass minpass]
            yd = [above above NaN below below];
        else
        % 'ydata': [ maxpass maxpass NaN minpass minpass NaN ...
        %                    maxpass maxpass NaN minpass minpass])
            yd = [above above NaN below below ...
                   NaN above above NaN below below ];
        end
        L1.ydata = yd;
        Rp.value = newRp;
        ylim = [below above];
        dyl = (ylim(2)-ylim(1))*.15;
        ax.ylimPassband = ylim + [-dyl/2 dyl/2];
    else
        Rp.value = newRp;
        RChangeMeas(1,type,Fs,Lresp,L1,L2,f1,f2,f3,f4,Rp,Rs)
    end
 
%---------------------------------------------------------------------
% fdbutter('Rsdrag',type,ind)
%  drag callback of L2 - stopband line
%  Inputs:
%     type - band configuration 1==low, 2=high, 3=pass, 4=stop
%     ind - index of segment being dragged
case 'Rsdrag'
    type = varargin{2};
    ind = varargin{3};
  
    yd = L2.ydata;
    newRs = -yd(ind);
    if newRs < 0
        newRs = 0;
    end
    if minOrdCheckbox.value
        switch type
        case {1,2,4}
            L2.ydata = [-newRs -newRs];
        case 3
            L2.ydata = [-newRs -newRs NaN -newRs -newRs];
        end
        set(Rs,'value', newRs)
    else
        set(Rs,'value', newRs)
        RChangeMeas(0,type,Fs,Lresp,L1,L2,f1,f2,f3,f4,Rp,Rs)
    end
    
end  % of function switch-yard

%---------------------------------------------------------------------
% -------- LOCAL FUNCTIONS START HERE  --------- 
%---------------------------------------------------------------------

function sethelp
global minOrdCheckbox bandpop 
global passframe stopframe passframe1 stopframe1
global f1 f2 f3 f4 w3db1 w3db2 order Rp Rm pbspecs sbspecs
global Rpm w3db1m w3db2m
global ax L1 L2 Lresp
global Fs

% disp('setting help ... stub')

function [n,w3db] = estimateOrder(type,Rp,Rs,Fs,f1,f2,f3,f4)
% [n,w3db] = estimateOrder(type,Rp,Rs,Fs,f1,f2,f3,f4)
%   estimate filter order
%   takes the specifications as given by the input
%   parameters and estimates the order and 3db frequencies
%   needed to meet those specifications.
% Inputs:
%  type - 1,2,3,4 specifies band configuration
%  Rp, Rs passband, stopband ripple
%  Fs - sampling frequency
%  f1,f2  first two frequencies in ascending order
%  f3,f4 only needed if type == 3 or 4, remaining frequencies
%    f1,f2,f3,f4 are assumed between 0 and Fs/2 on input
% Outputs:
%  n - filter order
%  w3db - filter band edges for butter, normalized to range [0...1]

if type == 1    % low pass
    Wp = f1;   Ws = f2;
elseif type == 2   % high pass
    Wp = f2;   Ws = f1;
elseif type == 3   % band pass
    Wp = [f2 f3];   Ws = [f1 f4];
elseif type == 4   % band stop
    Wp = [f1 f4];   Ws = [f2 f3];
end
[n,w3db] = buttord(Wp*2/Fs,Ws*2/Fs,Rp,Rs);
w3db = w3db(:)';   % make it a row 
 

function yd = passbandlimits(Rp)
% return ydata = [minpass maxpass] of passband 
%   given Rp decibels of ripple in passband (with maximum 1 in linear scale)
    above = 0; below = -Rp;
    yd = [below above];


function [maxpass,minpass,minstop] = getMagMeasurements(ff,Hlog,type,...
                                          f1,f2,f3,f4);
%getMagMeasurements
% Finds passband and stopband ripple for given band edges 
% given a filter's magnitude response
% Inputs:
%  ff - xdata of response
%  Hlog - magnitude of response at frequencies ff, in dB
%  type - band configuration (lp = 1, hp = 2, bp = 3, bs = 4)
%  f1, f2, f3, f4 - band edges (f3 and f4 ignored if type < 3)
%           in same units as ff
% Output:
%  fm - 2 or 4 element frequency vector in ascending order
    switch type
    case 1   % lowpass
        passInd = find(ff<=f1);
        stopInd = find(ff>=f2);
    case 2   % highpass
        stopInd = find(ff<=f1);
        passInd = find(ff>=f2);
    case 3   % bandpass
        stopInd = find((ff<=f1)|(ff>=f4));
        passInd = find((ff>=f2)&(ff<=f3));
    case 4   % bandstop
        passInd = find((ff<=f1)|(ff>=f4));
        stopInd = find((ff>=f2)&(ff<=f3));
    end
    maxpass = max(Hlog(passInd));
    minpass = min(Hlog(passInd));
    minstop = max(Hlog(stopInd));

function [fm,Rs] = getFreqMeasurements(ff,Hlog,type,Rp,Rs)
%getFreqMeasurements
% Finds band edges for a given passband and stopband ripple
% given a filter's magnitude response
% Inputs:
%  ff - xdata of response (assumed a column vector)
%  Hlog - magnitude of response at frequencies ff, in dB
%               (assumed a column vector)
%  type - band configuration (lp = 1, hp = 2, bp = 3, bs = 4)
%  Rp - passband ripple, in dB
%  Rs - stopband attenuation in dB
% Output:
%  fm - 2 or 4 element frequency vector in ascending order
%  Rs - changed value of Rs in case given Rs is unattainable
    passInd = find(Hlog(:)>=-Rp);
    stopInd = find((Hlog(:)<=-Rs) & (Hlog(:)<-Rp));
    switch type
    case 1   % lowpass
        if isempty(stopInd)
           [Rs1,stopInd] = min(Hlog);  Rs1 = -Rs1;
           if Rs1<Rs
               Rs = Rs1;
           else
               stopInd = passInd(end)+1;
               Rs = -Hlog(stopInd);
           end
        end
        fm = ff([passInd(end); stopInd(1)]);
    case 2   % highpass
        if isempty(stopInd)
           [Rs1,stopInd] = min(Hlog);  Rs1 = -Rs1;
           if Rs1<Rs
               Rs = Rs1;
           else
               stopInd = passInd(1)-1;
               Rs = -Hlog(stopInd);
           end
        end
        fm = ff([stopInd(end); passInd(1)]);
    case 3   % bandpass
        if isempty(stopInd)
           [Rs1,stopInd1] = min(Hlog(1:passInd(1)));    Rs1 = -Rs1;
           [Rs2,stopInd2] = min(Hlog(passInd(end):end));  Rs2 = -Rs2;
           Rs1 = max([Rs1 Rs2]);
           if Rs1<Rs
               Rs = Rs1;
               stopInd = [stopInd1 stopInd2];
           else
               stopInd = [passInd(1)-1 passInd(end)+1];
               Rs = -max(Hlog(stopInd));
           end
        end
        %  index of beginning of upper stop band
        i = min(find(stopInd>=passInd(1))); 
        if isempty(i)
            i = length(stopInd);
            if i == 1
                stopInd = [1 2];
                i = 2;
            end
        end
        fm = ff([stopInd(i-1); passInd([1 end]); stopInd(i)]);
    case 4   % bandstop
        if isempty(stopInd)
           [Rs1,stopInd] = min(Hlog);    Rs1 = -Rs1;
           if Rs1<Rs
               Rs = Rs1;
           else
               i = min(find(passInd>=stopInd(1)));
               stopInd = [passInd(i-1)+1 passInd(i)-1];
               Rs = -max(Hlog(stopInd))
           end
        end
        %  index of beginning of upper pass band
        i = min(find(passInd>=stopInd(1)));
        fm = ff([passInd(i-1); stopInd(1); stopInd(end); passInd(i)]);
    end

function fChangeMeas(i,type,Fs,Lresp,L1,L2,f1,f2,f3,f4,Rp,Rs)
%fChangeMeas - interactively track response when a frequency has changed
%  This function is meant to be called when the user
%     a) changes the f1,f2,f3, or f4 measurement
%  or b) drags the f1,f2,f3, f4 value
%  Xdata, Ydata, values of other frequencies and Rp and Rs are updated as
%  needed.
%  Inputs:
%    i - which frequency has just changed
%    type - band configuration (1 = lp, 2 = hp, 3 = bp, 4 = bs)
%    Fs - sampling frequency
%    Lresp, L1, L2 - fdline objects
%    f1,f2,f3,f4,Rp,Rs - fdmeas objects

    ff = Lresp.xdata;
    delf = ff(2)-ff(1);  % frequency spacing
    Hlog = Lresp.ydata;
    switch type
    case 1 % lowpass
        if i == 1     % passband edge
            ind = round(f1.value/delf) + 1;
            f1.value = ff(ind);
            Rp.value = -Hlog(ind);
            passbandChange = 1;
        else  % stopband edge
            ind = round(f2.value/delf) + 1;
            f2.value = ff(ind);
            set(Rs,'value', -Hlog(ind))
            passbandChange = 0;
        end
    case 2  % highpass
        if i == 1     % stopband edge
            ind = round(f1.value/delf) + 1;
            f1.value = ff(ind);
            set(Rs,'value', -Hlog(ind))
            passbandChange = 0;
        else  % passband edge
            ind = round(f2.value/delf) + 1;
            f2.value = ff(ind);
            Rp.value = -Hlog(ind);
            passbandChange = 1;
        end
    case 3 % bandpass
        switch i
        case 1     % lower stopband edge
            ind = round(f1.value/delf) + 1;
            f1.value = ff(ind);
            set(Rs,'value', -Hlog(ind))
            passbandChange = 0;
        case 2     % lower passband edge
            ind = round(f2.value/delf) + 1;
            f2.value = ff(ind);
            Rp.value = -Hlog(ind);
            passbandChange = 1;
        case 3     % upper passband edge
            ind = round(f3.value/delf) + 1;
            f3.value = ff(ind);
            Rp.value = -Hlog(ind);
            passbandChange = 1;
        case 4     % upper stopband edge
            ind = round(f4.value/delf) + 1;
            f4.value = ff(ind);
            set(Rs,'value', -Hlog(ind))
            passbandChange = 0;
        end
    case 4 % bandstop
        switch i
        case 1     % lower passband edge
            ind = round(f1.value/delf) + 1;
            f1.value = ff(ind);
            Rp.value = -Hlog(ind);
            passbandChange = 1;
        case 2     % lower stopband edge
            ind = round(f2.value/delf) + 1;
            f2.value = ff(ind);
            set(Rs,'value', -Hlog(ind))
            passbandChange = 0;
        case 3     % upper stopband edge
            ind = round(f3.value/delf) + 1;
            f3.value = ff(ind);
            set(Rs,'value', -Hlog(ind))
            passbandChange = 0;
        case 4     % upper passband edge
            ind = round(f4.value/delf) + 1;
            f4.value = ff(ind);
            Rp.value = -Hlog(ind);
            passbandChange = 1;
        end
    end
    if type >= 3  % find new freq. values in bp, bs case
        fm = getFreqMeasurements(ff,Hlog,type,Rp.value,Rs.value);
        switch type
        case 3 % bandpass - update f2 f3
            if passbandChange
                f2.value = fm(2);
                f3.value = fm(3);
            else
                f1.value = fm(1);
                f4.value = fm(4);
            end
        case 4 % bandstop - update f1 f4
            if passbandChange
                f1.value = fm(1);
                f4.value = fm(4);
            else
                f2.value = fm(2);
                f3.value = fm(3);