fdfirls.m

matlabDigitalSigalProcess内有文件若干

function varargout = fdfirls(varargin)
%fdfirls  Firls - Least Squares Module for filtdes.

%   Author: T. Krauss
%   Copyright (c) 1988-98 by The MathWorks, Inc.
%   $Revision: 1.5 $

% Change this global to static when available:
%  following are common to several modules
global minOrdCheckbox bandpop order pbspecs sbspecs pbmeasures sbmeasures
global passframe stopframe passframe1 stopframe1
%  following static to fdfirls
global f1 f2 f3 f4 wt1 wt2 Rp Rs 
global Rpm Rsm wt1m wt2m
global order0 order1
global ax L1 L2 Lresp
global Fs

switch varargin{1}
case 'init'

    filt = varargin{2};
    Fs = filt.Fs;

    [setOrderFlag_init, type_init, f_init, Rp_init, Rs_init, ...
      order_init, wt_init] = initSpecs(filt);
        
    [minOrdCheckbox bandpop order pbspecs sbspecs ...
         pbmeasures sbmeasures passframe stopframe ...
         passframe1 stopframe1 ax Lresp L1 L2 order1 ...
         L3_1 L3_2] = fdutil('commonObjects');
    order0 = order;
    order = order1;
        
    if ~strcmp(bandpop.userdata,'fdfirls')
       % this module is not current
        bandpop.userdata = 'fdfirls';
            
        minOrdCheckbox.callback = 'fdfirls(''checkbox'')';
        set(minOrdCheckbox,'visible','on','value',0)
        
        bandpop.callback = 'fdfirls(''newtype'')';         
        set(pbspecs(1),'callback','fdfirls(''fchange'')','range',[0 Fs/2])
        set(pbspecs(2),'callback','fdfirls(''fchange'')','range',[0 Fs/2])
        set(sbspecs(1),'callback','fdfirls(''fchange'')','range',[0 Fs/2])
        set(sbspecs(2),'callback','fdfirls(''fchange'')','range',[0 Fs/2])

        f3 = pbspecs(2);
        f4 = sbspecs(2);
        Rp = pbspecs(3);
        Rs = sbspecs(3);
        
        Rpm = pbmeasures(1);
        set(Rpm,'label','Actual Rp','format','%1.4g','visible','on')
        Rsm = sbmeasures(1);
        set(Rsm,'label','Actual Rs','format','%1.4g','visible','on')
        wt1m = pbmeasures(2);
        wt2m = sbmeasures(2);
        set(wt1m,'label','Weight','format','%1.4g')
        set(wt2m,'label','Weight','format','%1.4g')
        
        fdutil('changeToText',[pbmeasures(1:2); sbmeasures(1:2)])

        set(ax,'title','Frequency Response','xlimbound',[0 Fs/2],...
              'xlim',[0 Fs/2],...
              'xlabel','Frequency',...
              'ylabel','Magnitude (dB)');
              
        sethelp
            
        co = get(0,'defaultaxescolororder');
        lo = get(0,'defaultaxeslinestyleorder');
        
        % response line needs to be at bottom of stacking order:
        set(Lresp,'buttondownfcn','fdfirls(''LrespDown'')')
                       
        set(L1,'buttondownfcn','fdfirls(''L1down'')',...
               'buttonupfcn','fdfirls(''L1up'')');
    
        set(L2,'buttondownfcn','fdfirls(''L2down'')',...
               'buttonupfcn','fdfirls(''L2up'')');
    end    
    set(minOrdCheckbox,'visible','on')
    set(Rpm,'visible','on')
    set(Rsm,'visible','on')
    pbspecs(1).visible = 'on';
    pbspecs(3).visible = 'on';
    sbspecs(1).visible = 'on';
    sbspecs(3).visible = 'on';
    set(pbmeasures(3),'visible','off')
    set(sbmeasures(3),'visible','off')
    set(L3_1,'visible','off')
    set(L3_2,'visible','off')
    
    fdfirls('newfilt',setOrderFlag_init,type_init,f_init,Rp_init,...
               Rs_init,wt_init,order_init)
    
    [filt, errstr] = fdutil('callModuleApply',...
                               'fdfirls',filt,'');
    varargout{1} = filt;
    varargout{2} = errstr;

case 'apply'
    filt = varargin{2};
    msg = varargin{3};
    if strcmp(msg,'motion')
        Lresp.erasemode = 'xor';
    end
    
    setOrderFlag = 1;
    type = bandpop.value;
    
    % setup measurements:
    set(order1,'visible','off')
    set(wt1m,'visible','off')
    set(wt2m,'visible','off')
    
    % DESIGN FILTER!!!!
    n = order.value;
    switch type
    case 1   % lowpass
        f = [0 f1.value f2.value Fs/2]*2/Fs;
        m = [1 1 0 0];
        wt = [wt1.value wt2.value];
    case 2
        f = [0 f1.value f2.value Fs/2]*2/Fs;
        m = [0 0 1 1];
        wt = [wt2.value wt1.value];
    case 3
        f = [0 f1.value f2.value f3.value f4.value Fs/2]*2/Fs;
        m = [0 0 1 1 0 0];
        wt = [wt2.value wt1.value wt2.value];
    case 4
        f = [0 f1.value f2.value f3.value f4.value Fs/2]*2/Fs;
        m = [1 1 0 0 1 1];
        wt = [wt1.value wt2.value wt1.value];
    end
    if n>1000
        [continueFlag,errstr] = fdutil('largeWarning',n,msg);
        if ~continueFlag
            varargout{1} = filt;
            varargout{2} = errstr;
            return
        end
    end
    if ((type==2)|(type==4)) & rem(n,2)
        error('For highpass and bandstop filters, order must be even.')
    end
    b = firls(order.value,f,m,wt);
    a = 1;
    
    nfft = filtdes('nfft');
    [H,ff] = freqz(b,a,nfft,Fs);
    Hlog = 20*log10(abs(H));
    set(Lresp,'xdata',ff,'ydata',Hlog);

    [maxpass,minpass,minstop] = getMeasurements(ff,Hlog,bandpop.value,...
                                          f1.value,f2.value,f3.value,f4.value);
    Rpm.value = maxpass - minpass;
    Rsm.value = -minstop;
    ylim = [minpass maxpass];
    dyl = (ylim(2)-ylim(1))*.15;
    ax.ylimPassband = ylim + [-dyl/2 dyl/2];

    specstruc.setOrderFlag = 1;
    specstruc.type = bandpop.value;
    specstruc.f = f;
    specstruc.m = m;
    specstruc.Rp = Rpm.value;
    specstruc.Rs = Rsm.value;
    specstruc.wt = wt;
    specstruc.order = order.value;
    
    filt.tf.num = b;
    filt.tf.den = a;
    filt.specs.fdfirls = specstruc;
    filt.zpk = [];
    filt.ss = [];
    filt.sos = [];

    varargout{1} = filt;
    varargout{2} = '';
    
    switch type
    case 1   % lowpass
        set(L1,'xdata',[f(1:2) NaN f(1:2)]*Fs/2,...
               'ydata',[maxpass maxpass NaN minpass minpass])
        set(L2,'xdata',[f(3:4)]*Fs/2,'ydata',[minstop minstop])
    case 2   % highpass
        set(L1,'xdata',[f(3:4) NaN f(3:4)]*Fs/2,...
               'ydata',[maxpass maxpass NaN minpass minpass ])
        set(L2,'xdata',[f(1:2)]*Fs/2,'ydata',[ minstop minstop])
    case 3   % bandpass
        set(L1,'xdata',[f(3:4) NaN f(3:4)]*Fs/2,...
               'ydata',[maxpass maxpass NaN minpass minpass])
        set(L2,'xdata',[f(1:2) NaN f(5:6)]*Fs/2,...
               'ydata',[ minstop minstop NaN minstop minstop])
    case 4   % bandstop
        set(L1,'xdata',[f(1:2) NaN f(1:2) NaN f(5:6) NaN f(5:6)]*Fs/2,...
               'ydata',[ maxpass maxpass NaN minpass minpass NaN ...
                         maxpass maxpass NaN minpass minpass])
        set(L2,'xdata',[f(3:4)]*Fs/2,'ydata',[minstop minstop])
    end
    
    for i=1:3
        set(pbmeasures(i),'enable','on')
        set(sbmeasures(i),'enable','on')
    end
    set(order1,'enable','on')

    if strcmp(msg,'motion')
        return
    end
    
    Lresp.linestyle = '-';
    L1.linestyle = '-';
    L2.linestyle = '-';

case 'revert'
    filt = filtdes('filt');
    oldtype = filt.specs.fdfirls.type;
    % need to restore filter type
    setOrderFlag = filt.specs.fdfirls.setOrderFlag;
    f = filt.specs.fdfirls.f;
    Rpass = filt.specs.fdfirls.Rp;
    Rstop = filt.specs.fdfirls.Rs;
    wt = filt.specs.fdfirls.wt;
    N = filt.specs.fdfirls.order;
    fdfirls('newfilt',setOrderFlag,oldtype,f,Rpass,Rstop,wt,N)
    
    Rpm.value = L1.ydata(1)-L1.ydata(4);
    Rsm.value = -L2.ydata(1);
    
    Lresp.linestyle = '-';
    L1.linestyle = '-';
    L2.linestyle = '-';
    for i=1:3
        set(pbmeasures(i),'enable','on')
        set(sbmeasures(i),'enable','on')
    end
    set(order1,'enable','on')
    
case 'help'
    str = fdhelpstr('fdfirls');
    varargout{1} = str{2};

case 'description'
    varargout{1} = 'Least Squares FIR';

case 'Fs'
% Sampling frequency has changed

%   The filter spec does not depend on the sampling frequency
    filt = varargin{2};
    varargout{1} = filt;
    
% update various lines and specifications:
    oldFs = varargin{3};
    Fs = filt.Fs;  % new Fs
    f1.value = f1.value*Fs/oldFs;
    f2.value = f2.value*Fs/oldFs;
    f3.value = f3.value*Fs/oldFs;
    f4.value = f4.value*Fs/oldFs;
    f1.range = f1.range*Fs/oldFs;
    f2.range = f2.range*Fs/oldFs;
    f3.range = f3.range*Fs/oldFs;
    f4.range = f4.range*Fs/oldFs;

    ax.xlimbound = [0 Fs/2];
    ax.xlim = ax.xlim*Fs/oldFs;
    ax.xlimpassband = ax.xlimpassband*Fs/oldFs;
    L1.xdata = L1.xdata*Fs/oldFs;
    L2.xdata = L2.xdata*Fs/oldFs;
    Lresp.xdata = Lresp.xdata*Fs/oldFs;
    
%---------------------------------------------------------------------
% -------- following cases are module specific --------- 
%---------------------------------------------------------------------

case 'dirty'
% fdfirls('dirty')
% Callback of a spec ... to change appearance so user can
%  see that the specs and the currently designed filter
%  don't match
    if Lresp.linestyle ~= ':'
        Lresp.linestyle = ':';
        L1.linestyle = ':';
        L2.linestyle = ':';
    end
    for i=1:3
        set(pbmeasures(i),'enable','off')
        set(sbmeasures(i),'enable','off')
    end
    set(order1,'enable','off')

%---------------------------------------------------------------------
% fdfirls('fchange')
%  Callback when a frequency has changed
%  Need to update ranges and lines L1, L2
case 'fchange'
    xd1 = L1.xdata;
    xd2 = L2.xdata;
    f = [f1.value,f2.value,f3.value,f4.value];
    % update lines and find i, index of frequency which changed
    %  (i = 1, 2, 3, 4 for f1, f2, f3, or f4)
    i = [];
    type = bandpop.value;
    switch type
    case 1
        if xd1(2) ~= f(1)
            xd1([2 5]) = f(1);
            i = 1;
        else
            xd2(1) = f(2);
            i = 2;
        end
    case 2
        if xd1(2) ~= f(1)
            xd1(2) = f(1);
            i = 1;
        else
            xd2([1 4]) = f(2);
            i = 2;
        end
    case 3  % bandpass
      % L1 xdata = [f2 f3 NaN f2 f3]
      % L2 xdata = [0 f1 NaN f4 Fs/2]
        if xd2(2) ~= f(1)
            xd2(2) = f(1);
            i = 1;
        elseif xd1(1) ~= f(2)
            xd1([1 4]) = f(2);
            i = 2;
        elseif xd1(2) ~= f(3)
            xd1([2 5]) = f(3);
            i = 3;
        elseif xd2(4) ~= f(4)
            xd2(4) = f(4);
            i = 4;
        end
    case 4
      % L1 xdata = [0 f1 NaN 0 f1 NaN f4 Fs/2 NaN f4 Fs/2]
      % L2 xdata = [f2 f3]
        if xd1(2) ~= f(1)
            xd1([2 5]) = f(1);
            i = 1;
        elseif xd2(1) ~= f(2)
            xd2(1) = f(2);
            i = 2;
        elseif xd2(2) ~= f(3)
            xd2(2) = f(3);
            i = 3;
        elseif xd1(7) ~= f(4)
            xd1([7 10]) = f(4);
            i = 4;
        end
    end
    [xd1,xd2] = fdutil('validateBands',xd1,xd2,...
                       type,i,f,f1,f2,f3,f4,Fs);

    dirty = 0; % Flag used to indicate that specs have changed.
    
    % If there's a NaN in the data being compared the if-statement will
    % fail (return false) even tough the two variables are equal.  So,
    % remove the NaNs before comparing.
    nonNaNindx_L1 = ~isnan(L1.xdata);
    nonNaNindx_xd1 = ~isnan(xd1);
    if ~isequal(xd1(nonNaNindx_xd1),L1.xdata(nonNaNindx_L1)) 
        % only set if changed
        L1.xdata = xd1;
        dirty = 1;
    end

    % If there's a NaN in the data being compared the if-statement will
    % fail (return false) even tough the two variables are equal.  So,
    % remove the NaNs before comparing.
    nonNaNindx_L2 = ~isnan(L2.xdata);
    nonNaNindx_xd2 = ~isnan(xd2);
    if ~isequal(xd2(nonNaNindx_xd2),L2.xdata(nonNaNindx_L2)) 
        % only set if changed
        L2.xdata = xd2;
        dirty = 1;
    end

    if dirty  % Specs changed; update passband limits; show dotted lines
        ax.xlimPassband = fdutil('xlimpassband',type,...
                         Fs,f(1),f(2),f(3),f(4));
        fdfirls('dirty')
    end
               
%---------------------------------------------------------------------
% fdfirls('Rpchange')
%  Callback when Rp has changed
%  Need to update line L1
case 'Rpchange'

    Rpass = Rp.value;
    dev = (10^(Rpass/20)-1)/(10^(Rpass/20)+1);
    above = 20*log10(1+dev); below = 20*log10(1-dev);

    type = bandpop.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;
    ylim = [below above];
    dyl = (ylim(2)-ylim(1))*.15;
    ax.ylimPassband = ylim + [-dyl/2 dyl/2];

    fdfirls('dirty')
    
%---------------------------------------------------------------------
% fdfirls('Rschange')
%  Callback when Rs has changed
%  Need to update line L2
case 'Rschange'

    Rstop = Rs.value;
    type = bandpop.value;
    if type ~= 3
        yd = [-Rstop -Rstop];
    else
        yd = [-Rstop -Rstop NaN -Rstop -Rstop];
    end
    L2.ydata = yd;

    fdfirls('dirty')

%---------------------------------------------------------------------
% fdremez('checkbox')
%  Callback of minimum order checkbox
case 'checkbox'
    
    set(minOrdCheckbox,'value',0)
    
    error('For FIRLS, you can''t select "Minimum Order".')
    
%---------------------------------------------------------------------
% fdfirls('newfilt',setOrderFlag,type,f,Rp,Rs,wt,order)
%  set values of specs objects  DOES NOT DESIGN FILTER
case 'newfilt'
    setOrderFlag = varargin{2};
    type = varargin{3};
    f = varargin{4};
    Rpass = varargin{5};
    Rstop = varargin{6};
    wt = varargin{7};
    N = varargin{8};

    co = get(0,'defaultaxescolororder');
    
    bandpop.value = type;  
    % save last value of bandpop in passframe userdata:
    passframe.userdata = type;

    set(order0,'visible','on')
    order = order0;