fdellip.m

有关matlab的电子书籍有一定的帮助希望有用

    ylim = [below above];
    dyl = (ylim(2)-ylim(1))*.15;
    ax.ylimPassband = ylim + [-dyl/2 dyl/2];
 
%---------------------------------------------------------------------
% fdellip('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
    switch type
    case {1,2,4}
        L2.ydata = [-newRs -newRs];
    case 3
        L2.ydata = [-newRs -newRs NaN -newRs -newRs];
    end
    set(Rs,'value', newRs)
    
end  % of function switch-yard

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

function sethelp
global minOrdCheckbox bandpop 
global passframe stopframe passframe1 stopframe1
global f1 f2 f3 f4 Fp1 Fp2 order Rp Rm pbspecs sbspecs
global Fs1m Fs2m
global ax L1 L2 Lresp
global Fs

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

function [n,Fpass] = estimateOrder(type,Rp,Rs,Fs,f1,f2,f3,f4)
% [n,Fpass] = 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 stopband edge 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
%  Fpass - filter pass band edges for ellip, 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,Fpass] = ellipord(Wp*2/Fs,Ws*2/Fs,Rp,Rs);
Fpass = Fpass(:)';   % 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 Fstop = getFreqMeasurements(ff,Hlog,type,Rp,Rs,Fpass)
%getFreqMeasurements
% Finds stop band edges for elliptic filter given passband and 
% stopband ripple and 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
%  Fpass - passband edges, normalized to [0 Fs/2]
% Output:
%  Fstop - 1 or 2 element frequency (column) vector in ascending order
    ff = ff(:);
    switch type
    case 1   % lowpass
        stopInd = find(Hlog(:)<=-Rs);
        Fstop = ff(stopInd(1));
    case 2   % highpass
        stopInd = find(Hlog(:)<=-Rs);
        Fstop = ff(stopInd(end));
        %passInd = find(Hlog(:)>-Rs);
        %Fstop = ff(passInd(1)-1);
    case 3   % bandpass
        passInd = find(Hlog(:)>=-Rp);
        stopInd1 = find(Hlog(passInd(1):-1:1)<=-Rs);
        stopInd2 = find(Hlog(passInd(end):end)<=-Rs);
        Fstop = ff([passInd(1)-stopInd1(1) passInd(end)+stopInd2(1)]);
    case 4   % bandstop
        stopInd = find(Hlog(:)<=-Rs);
        Fstop = ff(stopInd([1 end]));
    end
    
    
function updateLines(passbandChange,type,fchange_ind,f1,f2,f3,f4,Rp,Rs,Fs)
% assume values of f1,f2,f3,f4,Rp and Rs are correct now
% fchange_ind - vector of indices indicating which frequencies have
%               changed
    global L1 L2 ax
    
    f = getFrequencyValues(type,f1,f2,f3,f4,Fs)*Fs/2;
    
    if any(passbandChange==1)
        maxpass = 0;
        minpass = -Rp.value;
        % update L1 xdata and ydata
        switch type
        case 1   % lowpass
            set(L1,'xdata',[f(1:2) NaN f(1:2)],...
                   'ydata',[maxpass maxpass NaN minpass minpass])
        case 2   % highpass
            set(L1,'xdata',[f(3:4) NaN f(3:4)],...
                   'ydata',[maxpass maxpass NaN minpass minpass ])
        case 3   % bandpass
            set(L1,'xdata',[f(3:4) NaN f(3:4)],...
                   'ydata',[maxpass maxpass NaN minpass minpass])
        case 4   % bandstop
            set(L1,'xdata',[f(1:2) NaN f(1:2) NaN f(5:6) NaN f(5:6)],...
                   'ydata',[ maxpass maxpass NaN minpass minpass NaN ...
                             maxpass maxpass NaN minpass minpass])
        end
        ylim = [minpass maxpass];
        dyl = (ylim(2)-ylim(1))*.15;
        ax.ylimPassband = ylim + [-dyl/2 dyl/2];
        if length(f)==4
            f(6)=0;  % zeropad for call to xlimpassband
        end
        ax.xlimPassband = fdutil('xlimpassband',type,...
                             Fs,f(2),f(3),f(4),f(5));
    end
    if any(passbandChange==0)
        minstop = -Rs.value;
        % update L2 xdata and ydata
        switch type
        case 1   % lowpass
            set(L2,'xdata',[f(3:4)],'ydata',[minstop minstop])
        case 2   % highpass
            set(L2,'xdata',[f(1:2)],'ydata',[ minstop minstop])
        case 3   % bandpass
            set(L2,'xdata',[f(1:2) NaN f(5:6)],...
                   'ydata',[ minstop minstop NaN minstop minstop])
        case 4   % bandstop
            set(L2,'xdata',[f(3:4)],'ydata',[minstop minstop])
        end
    end
    %fdutil('updateRanges',fchange_ind,f1,f2,f3,f4)
    set(f1,'range',[0 Fs/2])
    set(f2,'range',[0 Fs/2])
    set(f3,'range',[0 Fs/2])
    set(f4,'range',[0 Fs/2])
    
function [f1,f2,f3,f4] = setupFrequencyObjects(type,module,pbobjects,...
                   sbobjects,f,Fs,ax,setValueFlag)

    if nargin<8
        setValueFlag = 1;
    end
    
    switch type
    case 1 % lowpass
        f1 = pbobjects(1);
        f2 = sbobjects(1);
        f3 = sbobjects(2);
        f4 = pbobjects(2);
        pbobjects(2).visible = 'off';
        sbobjects(2).visible = 'off';
        set(f1,'label','Fp')
        set(f2,'label','Fs')
    case 2 % highpass
        f1 = sbobjects(1);
        f2 = pbobjects(1);
        f3 = pbobjects(2);
        f4 = sbobjects(2);
        pbobjects(2).visible = 'off';
        sbobjects(2).visible = 'off';
        set(f1,'label','Fs')
        set(f2,'label','Fp')
    case 3 % bandpass
        f1 = sbobjects(1);
        f2 = pbobjects(1);
        f3 = pbobjects(2);
        f4 = sbobjects(2);
        pbobjects(2).visible = 'on';
        sbobjects(2).visible = 'on';
        set(f1,'label','Fs1')
        set(f2,'label','Fp1')
        set(f3,'label','Fp2')
        set(f4,'label','Fs2')
    case 4
        f1 = pbobjects(1);
        f2 = sbobjects(1);
        f3 = sbobjects(2);
        f4 = pbobjects(2);
        pbobjects(2).visible = 'on';
        sbobjects(2).visible = 'on';
        set(f1,'label','Fp1')
        set(f2,'label','Fs1')
        set(f3,'label','Fs2')
        set(f4,'label','Fp2')
    end
    if setValueFlag
        if type < 3
            set(f1,'value',f(2)*Fs/2,'range',[0 f(3)]*Fs/2)
            set(f2,'value',f(3)*Fs/2,'range',[f(2) 1]*Fs/2)
        else
            set(f1,'value',f(2)*Fs/2,'range',[0 f(3)*Fs/2])
            set(f2,'value',f(3)*Fs/2,'range',[f(2) f(4)]*Fs/2)
            set(f3,'value',f(4)*Fs/2,'range',[f(3) f(5)]*Fs/2)
            set(f4,'value',f(5)*Fs/2,'range',[f(4) 1]*Fs/2)
        end
        ax.xlimPassband = fdutil('xlimpassband',type,...
                             Fs,f1.value,f2.value,f3.value,f4.value);
    end
    set(f1,'range',[0 Fs/2])
    set(f2,'range',[0 Fs/2])
    set(f3,'range',[0 Fs/2])
    set(f4,'range',[0 Fs/2])
    
    f1.callback = [module '(''fchange'',1)'];
    f2.callback = [module '(''fchange'',2)'];
    f3.callback = [module '(''fchange'',3)'];
    f4.callback = [module '(''fchange'',4)'];
        

function setupMeasurementObjects(specStruc,Fstop,Fs,...
                  order1,sbmeasures,Fs1m,Fs2m)
%  set values of MEASUREMENTS objects ... assumes specStruct is current

    setOrderFlag = specStruc.setOrderFlag;
    type = specStruc.type;
    N = specStruc.order;
    
    set(Fs1m,'value',Fstop(1),'format','%1.4g')
    if type > 2
        set(Fs1m,'label','Actual Fs1')
        set(Fs2m,'visible','on','value',Fstop(2),...
            'format','%1.4g','label','Actual Fs2')
    else
        set(Fs1m,'label','Actual Fs')
        set(Fs2m,'visible','off')
    end
    
    if ~setOrderFlag  % estimate order
        set(order1,'visible','on')
        order1.value = N;
    else
        set(order1,'visible','off')
    end
    for i=1:2
        set(sbmeasures(i),'enable','on')
    end
    set(order1,'enable','on')


function specStruc = measureFilt(objSetupFlag,specStruc,Fs,order1,...
                  sbmeasures,Fs1m,Fs2m,ff,Hlog,L1,L2,ax)

    n = specStruc.order;    
    setOrderFlag = specStruc.setOrderFlag;
    Fpass = specStruc.Fpass*Fs/2;
    
    Fstop = getFreqMeasurements(ff,Hlog,specStruc.type,...
                             specStruc.Rp,specStruc.Rs,Fpass);
    if objSetupFlag
        setupMeasurementObjects(specStruc,Fstop,Fs,...
                  order1,sbmeasures,Fs1m,Fs2m)    
    else
        set(Fs1m,'value',Fstop(1))
        if specStruc.type > 2
            set(Fs2m,'value',Fstop(2))
        end
        if ~setOrderFlag
            order1.value = n;
        end
    end

    if setOrderFlag
        % Update stopband edges in specStruc frequency vector:
        switch specStruc.type
        case 1
            f = [0 Fpass Fstop Fs/2];
        case 2
            f = [0 Fstop Fpass Fs/2];
        case 3
            f = [0 Fstop(1) Fpass(:)' Fstop(2) Fs/2];
        case 4
            f = [0 Fpass(1) Fstop(:)' Fpass(2) Fs/2];
        end
        specStruc.f = f*2/Fs;
    else
        f = specStruc.f*Fs/2;
    end
    
    maxpass = 0;   minpass = -specStruc.Rp;   minstop = -specStruc.Rs;
    
    % update L1 and L2 xdata and ydata, pointers, dragmodes
    fdutil('setLines','fdellip',L1,L2,specStruc.setOrderFlag,...
                      specStruc.type,specStruc.f(:)',...
                       Fs,minpass,maxpass,minstop)

    % update [ylim/xlim] passband limits
    ylim = [minpass maxpass];
    dyl = (ylim(2)-ylim(1))*.15;
    set(ax,'ylimPassband',ylim + [-dyl/2 dyl/2]);
    ax.xlimPassband = fdutil('xlimpassband',specStruc.type,...
                              Fs,f(2),f(3),f(end-2),f(end-1));

function f = getFrequencyValues(type,f1,f2,f3,f4,Fs);            
if type < 3  % low or high pass
    f = [0 f1.value f2.value Fs/2]*2/Fs;
else
    f = [0 f1.value f2.value f3.value f4.value Fs/2]*2/Fs;
end

function [setOrderFlag_init, type_init, f_init, Rp_init, Rs_init, ...
         order_init, Fpass_init] = initSpecs(filt)
%initSpecs   Initial specifications for ellip filter, from
%            filt input
%Switches off of filt.currentModule and if it finds any of
% fdcheby1, fdbutter, fdcheby2, fdremez, fdfirls, or fdkaiser
% retains the type, order, band edge, and any other relevant
% parameters

    % first define default values
    setOrderFlag_init = 0;   % by default, estimate order 
    type_init = 1;  % 1=lowpass, 2=highpass, 3=bandpass, 4=bandstop
    f_init = [0 .2 .25 1];
    Rp_init = 3;
    Rs_init = 20;
    order_init = 30;
    Fpass_init = .22;
    switch filt.specs.currentModule
    case {'fdcheby1','fdbutter','fdcheby2','fdremez','fdkaiser','fdfirls'}
        s = eval(['filt.specs.' filt.specs.currentModule]);
        setOrderFlag_init = s.setOrderFlag;
        type_init = s.type;
        f_init = s.f;
        Rp_init = s.Rp;
        Rs_init = s.Rs;
        order_init = s.order;
        switch filt.specs.currentModule
        case 'fdcheby1'
            Fpass_init = s.Fpass;
        case 'fdbutter'
            Fpass_init = s.w3db;
        case 'fdcheby2'
            Fpass_init = s.Fstop;
        case {'fdremez','fdfirls'}
            switch s.type
            case {1,2}
                Fpass_init = s.f(3);
            case 3
                Fpass_init = s.f([2 5]);
            case 4
                Fpass_init = s.f(3:4);
            end
        case 'fdkaiser'
            Fpass_init = s.Wn;
        end
        if any(strcmp(filt.specs.currentModule,...
                 {'fdremez','fdkaiser','fdfirls'}))
            order_init = ceil(order_init/10); 
             %  FIR filters are much higher order than IIR
        end
    case 'fdellip'
        if isfield(filt.specs,'fdellip')
            setOrderFlag_init = filt.specs.fdellip.setOrderFlag;
            type_init = filt.specs.fdellip.type;
            f_init = filt.specs.fdellip.f;
            Rp_init = filt.specs.fdellip.Rp;
            Rs_init = filt.specs.fdellip.Rs;
            order_init = filt.specs.fdellip.order;
            Fpass_init = filt.specs.fdellip.Fpass;
        end
    end