filtdemo.m

matlabDigitalSigalProcess内有文件若干

% initialize variables
    freqs = get(FreqsHndl,'UserData');  
    Fs = freqs(1);                   % sampling frequency
    Wp = freqs(2);                   % passband frequency
    Ws = freqs(3);                   % stopband frequency
    ripples = get(RipplesHndl,'UserData');
    Rp = ripples(1);                 % passband ripple
    Rs = ripples(2);                 % stopband attenuation
    
% Estimate order !
    % 1 = REMEZ, 2 = FIRLS, 4 = BUTTER, 5 = CHEBY1, 6 = CHEBY2, 7 = ELLIP, 
    %  3 = KAISER
    method = get(methodHndl,'value');
    if method == 1
        % compute deviations and estimate order
        devs = [ (10^(Rp/20)-1)/(10^(Rp/20)+1)  10^(-Rs/20) ];
        n = remezord([Wp Ws],[1 0],devs,Fs);
        n = max(3,n);
        wn = -1;
        max_n = 2000;  reasonable_n = 500;
    elseif method == 3
        devs = [ (10^(Rp/20)-1)/(10^(Rp/20)+1)  10^(-Rs/20) ];
        alfa = -min(20*log10(devs));
        if alfa > 50,
            beta = .1102*(alfa - 8.7);
        elseif alfa >= 21,
            beta = .5842*((alfa-21).^(.4)) + .07886*(alfa-21);
        else
            beta = 0;
        end
        n = ceil((alfa - 8)/(2.285*(Ws-Wp)*2*pi/Fs));
        n = max(1,n);
        wn = beta;
        max_n = 2000;  reasonable_n = 500;
    elseif method == 4
        [n,wn]=buttord(Wp*2/Fs,Ws*2/Fs,Rp,Rs);
        n = max(1,n);
        max_n = 60;  reasonable_n = 30;
    elseif method == 5
        [n,wn]=cheb1ord(Wp*2/Fs,Ws*2/Fs,Rp,Rs);
        n = max(1,n);
        max_n = 30;  reasonable_n = 15;
    elseif method == 6
        [n,wn]=cheb2ord(Wp*2/Fs,Ws*2/Fs,Rp,Rs);
        n = max(1,n);
        max_n = 30;  reasonable_n = 15;
    elseif method == 7
        [n,wn]=ellipord(Wp*2/Fs,Ws*2/Fs,Rp,Rs);
        n = max(1,n);
        max_n = 25;  reasonable_n = 12;
    end

    % set order field
    if method == 2   % can't estimate !
        set(ord1Hndl,'Userdata',[-1 -1],'String','-')
    else
        set(ord1Hndl,'Userdata',[n wn],'String',num2str(n))
    end

% determine which order to use: estimated or specified
    estim = get(btn1Hndl,'UserData');
    order = get(ord1Hndl,'UserData'); % use estimated Filter order
    wn = order(2);    % use estimated cutoff even for "specified order" case
    order = order(1);
    if estim ~= 1
        order = get(ord2Hndl,'UserData'); % use specified Filter order
    elseif (n>max_n)|(isnan(n))|(isinf(n))
        set(ord2Hndl,'UserData',reasonable_n,'string',num2str(reasonable_n))
        set(btn1Hndl,'value',0,'userdata',2)  % switch radio buttons
        set(btn2Hndl,'value',1)
        order = get(ord2Hndl,'UserData'); % use specified Filter order
    end

% Create desired frequency response and weight vector for FIR case
    f = [0 Wp Ws Fs/2]/Fs*2;
    m = [1  1  0 0];
    devs = [ (10^(Rp/20)-1)/(10^(Rp/20)+1)  10^(-Rs/20) ];
    w = [1 1]*max(devs)./devs;

    drawnow

% Design filter
    if method==1,
      b = remez(order,f,m,w); a = 1;
      title_str = sprintf('Order %g FIR Filter designed with REMEZ',order);
    elseif method==2,
      b = firls(order,f,m,w); a = 1;
      title_str = sprintf('Order %g FIR Filter designed with FIRLS',order);
    elseif method==3, % Kaiser window
      b = fir1(order,(Wp+(Ws-Wp)/2)*2/Fs,kaiser(order+1,wn)); a = 1;
      title_str = sprintf('Order %g FIR Filter designed with FIR1 and KAISER',order);
    elseif method==4, % butterworth
      [b,a] = butter(order,wn);
      title_str = sprintf('Order %g Butterworth IIR Filter',order);
    elseif method==5, % chebyshev type I
      [b,a] = cheby1(order,Rp,wn);
      title_str = sprintf('Order %g Chebyshev Type I IIR Filter',order);
    elseif method==6, % chebyshev type II
      [b,a] = cheby2(order,Rs,wn);
      title_str = sprintf('Order %g Chebyshev Type II IIR Filter',order);
    elseif method==7, % elliptic
      [b,a] =ellip(order,Rp,Rs,wn);
      title_str = sprintf('Order %g Elliptic IIR Filter',order);
    end
    [H,F] = freqz(b,a,max( 2048,nextpow2(5*max(length(b),length(a))) ),Fs);
    H = 20*log10(abs(H));
    axes(freqzHndl)
    axhndlList = get(freqzHndl,'UserData');
    if isempty(axhndlList)   % first time - happens at initialization phase
        % initialize the axes
        above = 20*log10(1+devs(1));
        below = 20*log10(1-devs(1));
        hndl = plot(F,H,[0 Wp NaN 0 Wp],[above above NaN below below], ...
                    [Ws Fs/2],[-Rs -Rs]);
        set(gcf,'handlevisibility','callback')
        if length(a) > 1
            set(hndl(1),'color',colors(1,:),'userdata',[b; a])  % save filter 
            % coefficients in userdata of line
        else
            set(hndl(1),'color',colors(1,:),'userdata',b) 
        end
        set(hndl(2:3),'color',colors(min(size(colors,1),2),:),'linewidth',2)
        set(hndl(2),'buttondownfcn','filtdemo(''dragRp'',1)')
        set(hndl(3),'buttondownfcn','filtdemo(''dragRs'',1)')
        grid on
        axhndlList = hndl;
        set(freqzHndl,'UserData',axhndlList)
        set(freqzHndl,'xlim',[0 Fs/2],'ylim',[-Rs*(1.5) max(2*Rp,Rs*.1)])
        xlabel('Frequency (Hz)')
        ylabel('Magnitude (dB)')
    else
        if length(a) > 1
            set(axhndlList(1),'UserData',[b; a])
        else
            set(axhndlList(1),'UserData',[b])
        end 
        filtdemo('axesredraw')
    end
    title(title_str)

    set(gcf,'Pointer','arrow');
    return

elseif strcmp(action,'axesredraw'), % redraw axes
    set(gcf,'Pointer','watch');
    axHndl=gca;
    fhndlList=get(gcf,'Userdata');
    freqzHndl = fhndlList(1);
    FreqsHndl = fhndlList(4);
    RipplesHndl = fhndlList(5);
    viewHndl = fhndlList(7);
    btn1Hndl = fhndlList(8);
    btn2Hndl = fhndlList(9);
    ord1Hndl = fhndlList(10);
    ord2Hndl = fhndlList(11);
    methodHndl = fhndlList(14);

    colors = get(gca,'colororder'); 
    
    order = get(ord1Hndl,'UserData'); % estimated Filter order
    wn = order(2);
    order = order(1);
    if get(btn1Hndl,'userdata') ~= 1
        order = get(ord2Hndl,'UserData'); % use specified Filter order
    end
    freqs = get(FreqsHndl,'UserData');  
    Fs = freqs(1);                   % sampling frequency
    Wp = freqs(2);                   % passband frequency
    Ws = freqs(3);                   % stopband frequency
    ripples = get(RipplesHndl,'UserData');
    Rp = ripples(1);                 % passband ripple
    Rs = ripples(2);                 % stopband attenuation
    method = get(methodHndl,'value'); 

    % Create desired frequency response and weight vector
    f = [0 Wp Ws Fs/2]/Fs*2;
    m = [1  1  0 0];
    devs = [ (10^(Rp/20)-1)/(10^(Rp/20)+1)  10^(-Rs/20) ];
    w = [1 1]*max(devs)./devs;

    axes(freqzHndl)
    axhndlList = get(freqzHndl,'UserData');
    h = get(axhndlList(1),'userdata');
    if size(h,1) > 1
        b = h(1,:); a = h(2,:);
    else
        b = h; a = 1;
    end
    [H,F] = freqz(b,a,max( 2048,nextpow2(5*max(length(b),length(a))) ),Fs);
    H = 20*log10(abs(H));
    set(axhndlList(1),'Xdata',F,'Ydata',H);
    if method >= 1 & method <= 3      % FIR case
        above = 20*log10(1+devs(1)); below = 20*log10(1-devs(1));
    else
        above = 0; below = -Rp;
    end
    set(axhndlList(2),'Ydata',[above above NaN below below],...
                      'Xdata',[0 Wp NaN 0 Wp],'erasemode','normal')
    set(axhndlList(3),'Ydata',[-Rs -Rs],'XData',[Ws Fs/2],'erasemode','normal')
    viewmode = get(viewHndl,'value');  % 1 == full, 2 == passband, 3 == stopbnd 
    if viewmode == 1
        xlim = [0 Fs/2];
        ymin = min(22*log10(1-devs(1)),-Rs*(1.5));
        ymax = max(2*20*log10(1+devs(1)),Rs*.1);
    elseif viewmode == 2
        if method >= 1 & method <= 3      % FIR case
            ymax = max(max(H(find(F<Wp))),20*log10(1+devs(1)))*1.1;
            ymin = min(min(H(find(F<Wp))),20*log10(1-devs(1)))*1.1;
        else  % IIR
            ymax = max(max(H(find(F<Wp)))*1.1,.1*Rp);
            ymin = min(min(H(find(F<Wp)))*1.1,-1.1*Rp);
        end
        xlim = [0 Wp+.2*(Ws-Wp)];
    elseif viewmode == 3
        ymax = max(max(H(find(F>Ws)))*.9,-Rs*.9);
        ymin = min(max(H(find(F>Ws))), -Rs*1.6);
        xlim = [Ws-.2*(Ws-Wp) Fs/2];
    end
    set(freqzHndl,'xlim',xlim)
    if all(finite([ymin ymax]))
        set(freqzHndl,'ylim',[ymin ymax])
    else
        disp('filtdemo: Response is Inf or NaN')
        set(freqzHndl,'ylimmode','auto')
    end
    xlabel('Frequency (Hz)')
    ylabel('Magnitude (dB)')

    set(gcf,'Pointer','arrow');
    return

elseif strcmp(action,'info'),
    set(gcf,'pointer','arrow')
	ttlStr = get(gcf,'name');
	hlpStr1= [...
 ' This demo lets you design lowpass       '
 ' digital filters.  You can set the       '
 ' filter design method, the sampling      '
 ' frequency (Fsamp), the passband and     '
 ' stopband edge frequencies (Fpass and    '
 ' Fstop), the desired amount of ripple    '
 ' in the passband (Rpass) and attenua-    '
 ' tion in the stopband (Rs), and the      '
 ' filter order, by using the controls     '
 ' on the right of the figure.             '
 '                                         '
 ' If you want to enter in a filter order  '
 ' other than the one in the "Auto"        '
 ' field, enter a number in the "Set"      '
 ' field.  You can change between auto-    '
 ' matic order selection and setting your  '
 ' own filter order by clicking on the     '
 ' "Auto" and "Set" radio buttons.         '];
	hlpStr2 = [...
 ' The REMEZ, FIRLS and KAISER methods     '
 ' design linear phase Finite Impulse      '
 ' Response (FIR) filters.                 '
 '                                         '
 ' The BUTTER, CHEBY1, CHEBY2, and ELLIP   '
 ' methods design Infinite Impulse         '
 ' Response (IIR) filters.                 '
 '                                         '
 ' FIR filters generally require a much    '
 ' higher order than IIR. The "Auto"       '
 ' filter order option is not available    '
 ' for the FIRLS filter.                   '
 '                                         '
 ' Note that for IIR filters, the pass-    '
 ' band magnitude specification is         '
 ' between 0 and -Rp decibels, while for   '
 ' FIR filters, it is centered at          '
 ' magnitude 1 with equiripples.           '];

	hlpStr3 = [...
 ' The popup menu just above the "info"    '
 ' button lets you select a region of      '
 ' the frequency response to view.  You    '
 ' can zoom-in on the passband, the        '
 ' stopband, or look at both bands at      '
 ' once (the "full" view).                 '
 '                                         '
 ' You can also click and drag the Fpass   '
 ' and Fstop indicators on the figure to   '
 ' change Rpass, Rstop, Fpass or Fstop.    '
 '                                         '
 ' Filename: filtdemo.m                    '];

    myFig = gcf;

    helpfun(ttlStr,hlpStr1,hlpStr2,hlpStr3);
    return  % avoid fancy, self-modifying code which
    % is killing the callback to this window's close button
    % if you press the info button more than once.
    % Also, a bug on Windows MATLAB is killing the 
    % callback if you hit the info button even once!

    % Protect against gcf changing -- Change close button behind
    % helpfun's back
    ch = get(gcf,'ch');
    for i=1:length(ch),
      if strcmp(get(ch(i),'type'),'uicontrol'),
        if strcmp(lower(get(ch(i),'String')),'close'),
          callbackStr = [get(ch(i),'callback') ...
            '; filtdemo(''closehelp'',[' num2str(myFig) ' '  num2str(gcf) '])'];
          set(ch(i),'callback',callbackStr)
          return
        end
      end
    end
    return

elseif strcmp(action,'closehelp'),
    % Restore help's close button callback behind helpfun's back
    ch = get(s(2),'ch');
    for i=1:length(ch),
      if strcmp(get(ch(i),'type'),'uicontrol'),
        if strcmp(lower(get(ch(i),'String')),'close'),
          callbackStr = get(ch(i),'callback');
          k = findstr('; filtdemo(',callbackStr);
          callbackStr = callbackStr(1:k-1);
          set(ch(i),'callback',callbackStr)
          break;
        end
      end
    end
    ch = get(0,'ch');
    if ~isempty(find(ch==s(1))), figure(s(1)), end % Make sure figure exists

elseif strcmp(action,'getfilt')
    shh = get(0,'ShowHiddenHandles');
    set(0,'ShowHiddenHandles','on')    
    fhndlList=get(gcf,'Userdata');
    set(0,'ShowHiddenHandles',shh)
    freqzHndl = fhndlList(1);
    axhndlList = get(freqzHndl,'UserData');
    h = get(axhndlList(1),'userdata');
    if size(h,1) > 1
        num = h(1,:); den = h(2,:);
    else
        num = h; den = 1;
    end
else
  disp(sprintf( ...
     'filtdemo: action string ''%s'' not recognized, no action taken.',action))
end    % if strcmp(action, ...lose all