yhzfiltdes.m

信号实验常用的simulink模型和mfile,可直接在matlan下运行。

        
 %   else
        % make changes since filter spec has changed: dotted line, enable
        %  'design now' menu item
   %     set(ud.ht.magline,'linestyle',':')
   %     set(ud.ht.designmenu,'enable','on')
%    end

    setptr(fig,'arrow')
    set(fig,'windowbuttonmotionfcn',wbmf)  % restore motion function

case 'estimate'
% -------------------------------------------------------------------------
% [n,wn,max_n,reasonable_n] = yhzfiltdes('estimate',fig)
%   estimate filter order
%   this section takes the specifications as given by the current
%   figure status and estimates the order and special parameters
%   needed to meet those specifications.
%   It stores the resultant order and special parameters in
%   the userdata of the figure:
%      ud.specs.order.auto (n) is a scalar
%        for firls, n is -1 since you can't estimate the filter order
%      ud.specs.special.auto (wn) is a row containing the 'special' 
%      estimated parameters:
%        remez - weights   (one element for each band)
%        firls - weights   (one element for each band)
%        kaiser - beta of Kaiser window (scalar), followed by cutoff freq(s)
%        butter - 3dB frequency(ies)
%        cheby1 - passband edge(s)
%        cheby2 - stopband edge(s)
%        ellip  - passband edge(s)
%   Sets the static text 'string' to num2str(n)
%   Updates the special parameters dialog box, if it exists
%   If there are output arguments, return the order and the special
%   parameters:
%      [n,wn,max_n,reasonable_n] = yhzfiltdes('estimate');
%
    if nargin < 2
        fig = gcf;
    else
        fig = varargin{1};
    end

  % get handles
    ud=get(fig,'Userdata');

  % initialize variables

    FType = ud.specs.FType;    %   1 = lp, 2 = hp, 3 = bp, 4 = bs
    Fs = ud.specs.Fs;     % sampling frequency
    spec = ud.specs;
    fp = spec.f(1,:);
    fs = spec.f(2,:);
    %if FType > 2
    %    f3 = spec.f(3);
    %    f4 = spec.f(4);
    %end
    Rp = spec.Rp;
    Rs = spec.Rs;

% Estimate order !
    % 1 = REMEZ, 2 = FIRLS, 4 = BUTTER, 5 = CHEBY1, 6 = CHEBY2, 7 = ELLIP, 
    %  3 = KAISER
    firmethod = ud.specs.firmethod;
    iirmethod1 = ud.specs.iirmethod(1);
    iirmethod2 = ud.specs.iirmethod(2);
    ir = ud.specs.ir;
    f=ud.specs.f;
    if ir(1)  %method == 1
        % compute deviations and estimate order
        %if type == 1    % low pass
        %    dev = [ (10^(Rp/20)-1)/(10^(Rp/20)+1)  10^(-Rs/20) ];
        %    f = [f1 f2];  m = [1 0];
        %elseif type == 2   % high pass
        %    dev = [ 10^(-Rs/20) (10^(Rp/20)-1)/(10^(Rp/20)+1)];
        %    f = [f1 f2];  m = [0 1];
        %elseif type == 3   % band pass
        %    dev = [ 10^(-Rs/20) (10^(Rp/20)-1)/(10^(Rp/20)+1) 10^(-Rs/20)];
        %    f = [f1 f2 f3 f4];  m = [0 1 0];
        %elseif type == 4   % band stop
        %    dev = [(10^(Rp/20)-1)/(10^(Rp/20)+1) ...
        %               10^(-Rs/20) (10^(Rp/20)-1)/(10^(Rp/20)+1)];
        %    f = [f1 f2 f3 f4];  m = [1 0 1];
        %end
        %[n,fo,mo,wn] = remezord(f,m,dev,Fs);
        %max_n = 2000;  reasonable_n = 500;
        %wn = wn(:)';   % make it a row
        %if type == 2 | type == 4   % highpass or bandstop
        %    if rem(n,2)==1   % if n is odd, must increase
        %        n = n+1;
        %    end
        %end
        %if n<3, n=3; end   % n must be at least 3 for remez
    %elseif method == 2  % FIRLS
    %    if type == 1    % low pass
    %        dev = [ (10^(Rp/20)-1)/(10^(Rp/20)+1)  10^(-Rs/20) ];
    %    elseif type == 2   % high pass
    %        dev = [ 10^(-Rs/20) (10^(Rp/20)-1)/(10^(Rp/20)+1)];
    %    elseif type == 3   % band pass
    %        dev = [ 10^(-Rs/20) (10^(Rp/20)-1)/(10^(Rp/20)+1) 10^(-Rs/20)];
    %    elseif type == 4   % band stop
    %        dev = [(10^(Rp/20)-1)/(10^(Rp/20)+1) ...
    %                   10^(-Rs/20) (10^(Rp/20)-1)/(10^(Rp/20)+1)];
    %    end
    %    n = -1;   % can't estimate this baby
    %    wn = ones(size(dev))*max(dev)./dev;
   %     max_n = 2000;  reasonable_n = 500;
      if firmethod == 5   % KAISER
        dev = [ (10^(Rp/20)-1)/(10^(Rp/20)+1)  10^(-Rs/20) ];
        alfa = -min(20*log10(dev));
        if alfa > 50,
            beta = .1102*(alfa - 8.7);
        elseif alfa >= 21,
            beta = .5842*((alfa-21).^(.4)) + .07886*(alfa-21);
        else
            beta = 0;
        end
        if FType <= 2    % lowpass or highpass
            f1=min(f); f2=max(f);
            n = ceil((alfa - 8)/(2.285*(f2-f1)*2*pi/Fs));
            f = (f1+(f2-f1)/2);
         else  % bandpass 
            f1=f(2,1); f2=f(1,1); f3=f(1,2); f4=f(2,2);
            n = ceil((alfa - 8)/(2.285*min(f2-f1,f4-f3)*2*pi/Fs));
            f = [(f1+(f2-f1)/2) (f3+(f4-f3)/2)];
        end
        wn = [beta f];
        max_n = 2000;  reasonable_n = 500;
        %if FType == 2 | FType == 4   % highpass or bandstop
            if rem(n,2)==1   % if n is odd, must decrease
                n = n-1;
            end
         %end
      ud.specs.kaiserbeta=beta;
      ud.specs.firN=n+1;
      elseif firmethod ==6 %Frequency Sampling Method
         %???
      else 
         n=ud.specs.firN;
      end
     % ud.specs.firN=n;
   else %iir   
      if iirmethod1==1
      %elseif method == 4
        if FType == 1    % lowpass
           [n,wn]=buttord(fp*2/Fs,fs*2/Fs,Rp,Rs);
        elseif FType == 2   % highpass
           [n,wn]=buttord(fp*2/Fs,fs*2/Fs,Rp,Rs);
        elseif FType == 3   % bandpass
           [n,wn]=buttord(fp*2/Fs,fs*2/Fs,Rp,Rs);
    %    elseif type == 4   % bandstop
    %        [n,wn]=buttord([f1 f4]*2/Fs,[f2 f3]*2/Fs,Rp,Rs);
        end
       wn = wn*Fs/2;
        max_n = 60;  reasonable_n = 30;
      elseif iirmethod1 == 2
       % if type == 1    % lowpass
       %     [n,wn]=cheb1ord(f1*2/Fs,f2*2/Fs,Rp,Rs);
       % elseif type == 2   % highpass
       %     [n,wn]=cheb1ord(f2*2/Fs,f1*2/Fs,Rp,Rs);
       % elseif type == 3   % bandpass
       %     [n,wn]=cheb1ord([f2 f3]*2/Fs,[f1 f4]*2/Fs,Rp,Rs);
       % elseif type == 4   % bandstop
       %     [n,wn]=cheb1ord([f1 f4]*2/Fs,[f2 f3]*2/Fs,Rp,Rs);
       % end
       [n,wn]=cheb1ord(fp*2/Fs,fs*2/Fs,Rp,Rs);
        wn = wn*Fs/2;
        max_n = 30;  reasonable_n = 15;
      elseif iirmethod1 == 3
        %if type == 1    % lowpass
        %    [n,wn]=cheb2ord(f1*2/Fs,f2*2/Fs,Rp,Rs);
        %elseif type == 2   % highpass
        %    [n,wn]=cheb2ord(f2*2/Fs,f1*2/Fs,Rp,Rs);
        %elseif type == 3   % bandpass
        %    [n,wn]=cheb2ord([f2 f3]*2/Fs,[f1 f4]*2/Fs,Rp,Rs);
        %elseif type == 4   % bandstop
        %    [n,wn]=cheb2ord([f1 f4]*2/Fs,[f2 f3]*2/Fs,Rp,Rs);
        %end
        [n,wn]=cheb2ord(fp*2/Fs,fs*2/Fs,Rp,Rs);
        wn = wn*Fs/2;
        max_n = 30;  reasonable_n = 15;
    %elseif method == 7   % elliptic
    %    if type == 1    % lowpass
    %        [n,wn]=ellipord(f1*2/Fs,f2*2/Fs,Rp,Rs);
    %    elseif type == 2   % highpass
    %        [n,wn]=ellipord(f2*2/Fs,f1*2/Fs,Rp,Rs);
    %    elseif type == 3   % bandpass
    %        [n,wn]=ellipord([f2 f3]*2/Fs,[f1 f4]*2/Fs,Rp,Rs);
    %    elseif type == 4   % bandstop
    %        [n,wn]=ellipord([f1 f4]*2/Fs,[f2 f3]*2/Fs,Rp,Rs);
    %    end
    %    wn = wn*Fs/2;
    %    max_n = 25;  reasonable_n = 12;
       end

    % set order field
    ud.specs.iirorder = n;
    %ud.specs.special.auto = wn;
    %if method == 2   % can't estimate !
    %    set(ud.ht.ord1Hndl,'String','-')
    %else
    end
    if ir(2)
        set(ud.ht.iirorderHndl,'String',num2str(n))
     elseif firmethod==5
        label=sprintf('N:\nbeta:');
        set(ud.ht.firlabelHndl,'string',label); %'string','[N\nbeta]:');
        label=sprintf('%g\n%g',ud.specs.firN,beta);
        set(ud.ht.firNHndl,'string',label);%[num2str(n),' ',num2str(beta)]);
        %end
     else
        set(ud.ht.firlabelHndl,'string','N:');
        set(ud.ht.firNHndl,'string',num2str(ud.specs.firN));
    end   
    set(fig,'userdata',ud)

    %fdspdlg('updateparamdialog',fig)

    % output params?  Return estimates if so
    if nargout > 0
        varargout{1} = n;
    end
    if nargout > 1
        varargout{2} = wn;
    end
    if nargout > 2
        varargout{3} = max_n;
    end
    if nargout > 1
        varargout{4} = reasonable_n;
    end

case 'edit'
% -------------------------------------------------------------------------
% yhzfiltdes('edit')
%   callback of editbox for filter parameters (band edges and attenuations)
%
    ud = get(gcf,'UserData');
    FType = ud.specs.FType;    %%% 1 = lp, 2 = hp, 3 = bp, 4 = bs

    if FType<=2
       v = [ud.specs.f; ud.specs.Rp; ud.specs.Rs];   % old values
    else 
       v = [ud.specs.f(1,:)';ud.specs.f(2,:)'; ud.specs.Rp;ud.specs.Rs];
    end   
    str = get(gco,'String');
    ind = find(abs(str)<32);
    str(ind) = 32*ones(size(ind));
    str = [str 32*ones(size(str,1),1)];   % add a space after each number
    str = str';
    str = char(['[' str(:)' ']']);
    vv = eval(str,'-1')';     % new values
    if ~isempty(vv) & vv == -1
        vv = v;    % in case of error, restore old values
    end
    %if FType <= 2
    if isempty(vv) | any(imag(vv)) | any(vv<=0) %| length(vv)~=4  %| any(diff(vv(1:2))<=0) |             
       vv = v;
    end
    inpError = 0;
    if FType ==1 &  vv(1)>vv(2)
          inpError=1;
    end   
    if FType ==2 & vv(1)<vv(2)
          inpError =1;
    end
    if FType ==3 & ~(vv(3)<=vv(1) & vv(1)<=vv(2) & vv(2)<=vv(4))
          inpError =1;
    end
    if FType<3 & any(vv(1:2)>ud.specs.Fs/2)
       inpError=1;
    elseif FType==3 & any(vv(1:4)>ud.specs.Fs/2)
       inpError=1;
    end;      
    if inpError
       disp('输入错误, 请重新输入!');
       vv = v;
    end
    if FType<=2
        ud.specs.f = vv(1:2);
        ud.specs.Rp = vv(3);
        ud.specs.Rs = vv(4);
    else
       ud.specs.f=[vv(1:2)';vv(3:4)'];
       ud.specs.Rp = vv(5);
       ud.specs.Rs=vv(6);
    end   
       %if isempty(vv) | any(imag(vv)) | any(vv<=0) | ...
    %              any(diff(vv(1:4))<=0) | length(vv)~=6
    %        vv = v;
    %    end
    %    ud.specs.f = vv(1:4);
    %    ud.specs.Rp = vv(5);
    %    ud.specs.Rs = vv(6);
    %end
%    if spec.type <= 2
%        str = sprintf('%g\n%g\n\n%g\n%g',spec.f(1),spec.f(2),...
%          spec.Rp,spec.Rs);
%    else
%        str = sprintf('%g\n%g\n%g\n%g\n\n%g\n%g',spec.f(1),spec.f(2),...
%          spec.f(3),spec.f(4),spec.Rp,spec.Rs);
%    end
    set(gco,'String',Localfdspecstr(ud.specs))
    set(gcf,'userdata',ud)
    v = v(:);
    if any(vv~=v)
        yhzfiltdes('specchange')
    end
    
 case 'setFs'
% -------------------------------------------------------------------------
% yhzfiltdes('setFs',fig,dogrid)
%   callback to set sampling frequency of tool
%
    if nargin < 2
        fig = gcf;
    else
        fig = varargin{1};
    end

    if nargin < 3
        dogrid = 1;
    else
        dogrid = 0;
    end

    ud = get(fig,'userdata');

    %oldFs = ud.specs.Fs % old value of Fs sampling frequency
    oldFs = get(ud.ht.FsHndl,'userdata');  % old value of Fs sampling frequency
    str = get(ud.ht.FsHndl,'String');
    Fs = eval(str);
    if Fs<=0 | imag(Fs)~=0 | isempty(Fs)
       disp('输入错误, 请重新输入!');
       Fs=oldFs;
    end   
    %[Fs,err] = validarg(str,[0 Inf],[1 1],'sampling frequency');
    %Fs = ud.filt.Fs;     % new value (only different when the sampling frequency
    %                     %            is changed in the SPTool).
    %toolFs = get(ud.ht.FsHndl,'userdata');  % value currently in yhzfiltdes
    
    set(ud.ht.FsHndl,'userdata',Fs,'string',str);   % store new Fs
    ud.specs.Fs = Fs;
    if (Fs-oldFs)    % if Fs really changed,
        % update the tool here

        % update the band edges
        ud.specs.f = ud.specs.f*Fs/oldFs;
        
        % update the "Full View" limits
        ud.limits.xlim = ud.limits.xlim*Fs/oldFs;
        
        % update prefs string
        ud.prefs.Fs = str;

        if dogrid
          % make sure freq. grid spacing is proportionally the same
          ud.prefs.gridsize(1) = ud.prefs.gridsize(1)*Fs/oldFs;
          %if ~isempty(ud.tabfig)
          %   tabud = get(ud.tabfig,'userdata');
          %   if ~isempty(get(tabud.tabs(2),'userdata'))
          %      fdprefhand('populate',ud.tabfig,2,ud.prefs)
          %   end
          %end
        end
        
        % save the user data
        set(fig,'userdata',ud)

        % set string of filter specs edit box:
        set(ud.ht.specHndl,'string',Localfdspecstr(ud.specs)); 

        yhzfiltdes('drawspeclines',fig)

        set(ud.ht.ax1,'xlim',get(ud.ht.ax1,'xlim')*Fs/oldFs); 
        set(ud.ht.magline,...
              'xdata',get(ud.ht.magline,'xdata')*Fs/oldFs);  % magline

        figname = get(fig,'name');

        %if ~isempty(ud.tabfig)    % update settings figure
        %    % since Fs is on the first page, it has always been created
        %    fdprefhand('populate',ud.tabfig,1,ud.prefs)
        %end
        
        % update special parameters - both estimated and manually selected 
        firmethod = ud.specs.firmethod;
        iirmethod = ud.specs.iirmethod;
        %wn_est = ud.specs.special.auto;
        %wn_sel = ud.specs.special.manual;
        %if method == 3   % kaiser window
        %    wn_est(2:length(wn_est)) = wn_est(2:length(wn_est))*Fs/oldFs;
        %    wn_sel(2:length(wn_sel)) = wn_sel(2:length(wn_sel))*Fs/oldFs;
        %elseif method >= 4  % IIR filter