www.gusucode.com > signal 工具箱matlab源码程序 > signal/stepz.m

    function varargout = stepz(b,varargin)
%STEPZ  Step response of digital filter
%   [H,T] = STEPZ(B,A) computes the step response of the filter:
%
%               jw               -jw              -jmw
%        jw  B(e)    b(1) + b(2)e + .... + b(m+1)e
%     H(e) = ---- = ------------------------------------
%               jw               -jw              -jnw
%            A(e)    a(1) + a(2)e + .... + a(n+1)e
%
%   given numerator and denominator coefficients in vectors B and A.
%
%   [H,T] = STEPZ(SOS) computes the step response of the filter specified
%   using the second order sections matrix SOS. SOS is a Kx6 matrix, where
%   the number of sections, K, must be greater than or equal to 2. Each row
%   of SOS corresponds to the coefficients of a second order filter. From
%   the transfer function displayed above, the ith row of the SOS matrix
%   corresponds to [bi(1) bi(2) bi(3) ai(1) ai(2) ai(3)].
%
%   [H,T] = STEPZ(D) computes the step response of the digital filter, D.
%   You design a digital filter, D, by calling the <a href="matlab:help designfilt">designfilt</a> function.
%
%   In all cases the number of samples is chosen internally, and STEPZ
%   returns the response in column vector H and a vector of times (or
%   sample intervals) in T (T = [0 1 2 ...]').
%
%   [H,T] = STEPZ(...,N) computes the first N samples of the step response.
%   If N is a vector of integers, the step response is computed only at
%   those integer values (0 is the origin).
%
%   [H,T] = STEPZ(...,N,Fs) separates the time samples by T = 1/Fs seconds.
%   Fs is assumed to be in Hz.
%
%   STEPZ(...) with no output arguments plots the step response.
%
%   % Example 1:
%   %   Design a lowpass FIR filter with normalized cut-off frequency at
%   %   0.3 and show its step response.
%
%   b=fircls1(54,0.3,0.02,0.008);
%   stepz(b)
%
%   % Example 2:
%   %   Display the step response of a 3rd order Butterworth IIR filter.
%
%   [b,a] = butter(3,0.4);
%   stepz(b,a)
%
%   % Example 3:
%   %   Design a Butterworth highpass IIR filter, represent its coefficients
%   %   using second order sections, and display its step response.
%
%   [z,p,k] = butter(6,0.7,'high');
%   SOS = zp2sos(z,p,k);
%   stepz(SOS)
%
%   % Example 4:
%   %   Use the designfilt function to design a highpass IIR digital filter 
%   %   with order 8, passband frequency of 75 KHz, and a passband ripple 
%   %   of 0.2 dB. Sample rate is 200 KHz. Visualize the step response 
%   %   using 256 samples.
%  
%   D = designfilt('highpassiir', 'FilterOrder', 8, ...
%            'PassbandFrequency', 75e3, 'PassbandRipple', 0.2,...
%            'SampleRate', 200e3);
%
%   stepz(D,256)
%
%   See also IMPZ, FREQZ, ZPLANE, GRPDELAY, FVTOOL.

%   Copyright 1988-2013 The MathWorks, Inc.

narginchk(1,4)

isTF = true; % True if dealing with a transfer function

if all(size(b)>[1 1])
    % Input is a matrix, check if it is a valid SOS matrix
    if size(b,2) ~= 6
        error(message('signal:signalanalysisbase:invalidinputsosmatrix'));
    end
    isTF = false; % SOS instead of transfer function
    
    if nargin > 1
        N = varargin{1};
    else
        N = [];
    end
    if nargin > 2
        Fs = varargin{2};
    else
        Fs = 1;
    end
    % Cast to enforce precision rules
   N = signal.internal.sigcasttofloat(N,'double','stepz','N',...
     'allownumeric');
   Fs = signal.internal.sigcasttofloat(Fs,'double','stepz','Fs',...
     'allownumeric');
    
   % Checks if SOS is a valid numeric data input
   signal.internal.sigcheckfloattype(b,'','stepz','SOS');
end

% If input is a transfer function b,a -----
if isTF
    if nargin > 1
        a = varargin{1};
        if all(size(a)>[1 1])
            error(message('signal:signalanalysisbase:inputnotsupported'));
        end
    else
        a = 1;
    end
    
    if nargin > 2
        N = varargin{2};
    else
        N = [];
    end
    
    if nargin > 3
        Fs = varargin{3};
    else
        Fs = 1;
    end
    % Cast to enforce precision rules
    N = signal.internal.sigcasttofloat(N,'double','stepz','N',...
      'allownumeric');
    Fs = signal.internal.sigcasttofloat(Fs,'double','stepz','Fs',...
      'allownumeric');
      
    % Checks if B and A are valid numeric data inputs
    signal.internal.sigcheckfloattype(b,'','stepz','B');
    signal.internal.sigcheckfloattype(a,'','stepz','A');
end

% Compute time vector
M = 0;  NN = [];
if isempty(N)
    % if not specified, determine the length
    if isTF
        N = impzlength(b,a);
    else
        N  = impzlength(b);
    end
elseif length(N)>1  % vector of indices
    
    NN = round(N);
    N = max(NN)+1;
    M = min(min(NN),0);
end

t = (M:(N-1))'/Fs;

% Form input vector
x = ones(size(t));

if isTF
    s = filter(b,a,x);
else
    s = sosfilt(b,x);
end

if ~isempty(NN),
    s = s(NN-M+1);
    t = t(NN-M+1);
end

% Cast to enforce precision rules.
% Only cast if output is requested, otherwise, plot using double precision
% time vector. 
if nargout
    if isa(s,'single')
      t = single(t);
    end
    varargout = {s,t};
else
    timezplot(t,s,Fs,getString(message('signal:stepz:Step')));
end