显示MATLAB filter函数的源码

MATLAB 的 filter 函数的源码如下：

function [y,zf] = filter(b,a,x,zi,dim)
%FILTER 1-D digital filter.
%   Y = FILTER(B,A,X) filters the data in vector X with the filter described by
%   vectors A and B to create the filtered data Y.  The filter is a "Direct
%   Form II Transposed" implementation of the standard difference equation:
%
%     a(1)*y(n) = b(1)*x(n) + b(2)*x(n-1) + ... + b(nb+1)*x(n-nb)
%                         - a(2)*y(n-1) - ... - a(na+1)*y(n-na)
%
%   The initial states of the filter are assumed to be zero.
%
%   If A(1) is not equal to 1, FILTER normalizes the filter coefficients by A(1).
%
%   Y = FILTER(B,A,X,Z) gives access to initial and final conditions, Z.  If Z
%   is a scalar, FILTER applies the IC to both the input and output data.  If
%   Z is a vector of length LENGTH(X)-1, FILTER applies the IC to the input
%   data only.  If Z is the empty matrix [], FILTER applies no IC.
%
%   Y = FILTER(B,A,X,[],DIM) or Y = FILTER(B,A,X,Z,DIM) operates along the
%   dimension DIM. The DIM'th dimension of X must have length greater than 1.
%
%   [Y,ZF] = FILTER(B,A,X,Z,...) returns the final conditions, ZF. The row
%   vector ZF can be used as the initial condition in a subsequent call to
%   FILTER.
%
%   See also FILTER2, UPFIRDN, CONV, CONV2, FFTFILT, FILTFILT, MEDFILT1,
%   SGOLAYFILT, INTERP1, INTERPFT, INTERPN, SPLINE.
%   Filter Designer GUI:  FDESIGN.
%
%   Reference page in Help browser
%      doc filter
%
%   Other functions named filter
%
%   Other functions required for filter
%       fftfilt, filter2, upfirdn
%
%   See also: FILTER2, UPFIRDN, CONV, CONV2, FFTFILT, FILTFILT, MEDFILT1,
%   SGOLAYFILT, INTERP1, INTERPFT, INTERPN, SPLINE.

%   Copyright 1984-2019 The MathWorks, Inc.

% Check number of inputs
narginchk(3,5);

% Check order of inputs for backwards compatibility
if nargin==4 && isscalar(zi)
    dim = zi; zi = [];
end

% Check inputs
[rows,cols] = size(x);
if ~isnumeric(x)
    error(message('signal:filter:InvalidDataType'));
end
if ~isvector(a) || ~isvector(b)
    error(message('signal:filter:CoefficientsNotVectors'));
end
if (length(a) == 1) && (length(b) == 1)
    y = b * x;
    zf = [];
    return
end
if length(a) < 2 || length(b) < 1
    error(message('signal:filter:InvalidDimensionsOfCoefficients'));
end
if length(a) > length(b)
    b((end+1):length(a)) = 0;
elseif length(a) < length(b)
    a((end+1):length(b)) = 0;
end
if any(a==0)
    error(message('signal:filter:InvalidDenominator'));
end
if nargin < 5
    dim = find(size(x)>1,1);
    if isempty(dim), dim = 1; end
end
if isempty(zi)
    % Generate initial conditions
    zi = zeros(1,length(a)-1);
elseif isscalar(zi)
    % Expand scalar IC to vector
    zi = repmat(zi,1,length(a)-1);
else
    % Check IC dimensions
    if length(zi) ~= length(a)-1
        error(message('signal:filter:InvalidInitialConditions'));
    end
end

% Call filter implementation
if isequal(dim,1)
    [y,zf] = filter1D(b,a,x,zi);
elseif isequal(dim,2)
    [y,zf] = filter1D(b,a,x.',zi);
    y = y.';
else
    nDimsX = ndims(x);
    permOrder = [dim,1:dim-1,dim+1:nDimsX];
    xperm = permute(x,permOrder);
    sizperm = size(xperm);
    xperm = reshape(xperm,sizperm(1),[]);
    [y,zf] = filter1D(b,a,xperm,zi);
    y = reshape(y,sizperm);
    y = ipermute(y,permOrder);
    if ~isempty(zi)
        zf = reshape(zf,[length(a)-1 sizperm(2:end)]);
        zf = ipermute(zf,permOrder);
    end
end

%-------------------------------------------------------------------------
function [y,zf] = filter1D(b,a,x,zi)
%FILTER1D Perform 1-D filtering.

na = length(a);
nb = length(b);
nSamples = size(x,1);

% Delay line for filter
if na > 1
    z = zeros(nSamples,na-1);
    if ~isempty(zi)
        z(1,:) = zi(:).';
    end
else
    z = [];
end

% Allocate output
if isreal(x) && all(isfinite(x(:))) && all(isfinite(b(:))) && ...
   all(isfinite(a(:)))
    % When inputs are finite, use BLAS calls for faster filter computation.
    % Use real arithmetic when inputs and coefficients are real.
    y = zeros(size(x),class(x));
    if na > 1
        for k = 1:nSamples
            if nb > 1
                y(k) = b(1)*x(k) + sum(b(2:end).*z(k,:));
                z(k,:) = [y(k) z(k,1:end-1)];
            else
                y(k) = b*x(k);
            end
            if na > 1
                y(k) = y(k) - sum(a(2:end).*z(k,:));
            end
        end
    else
        y = filter(b,x);
    end
else
    % Use slower loop-based implementation
    y = filterScalarLoop(b,a,x,z);
end

% Compute final conditions
if nargout == 2 && ~isempty(z)
    zf = z(nSamples,:);
else
    zf = [];
end

%-------------------------------------------------------------------------
function y = filterScalarLoop(b,a,x,z)
%FILTERSCALARLOOP Perform filtering with scalar loop.

na = length(a);
nb = length(b);
nSamples = size(x,1);

% Allocate output
y = zeros(size(x),class(x));

% Filtering loop
for k = 1:nSamples
    if nb > 1
        y(k) = b(1)*x(k) + sum(b(2:end).*z(k,:));
        z(k,:) = [y(k) z(k,1:end-1)];
    else
        y(k) = b*x(k);
    end
    if na > 1
        y(k) = y(k) - sum(a(2:end).*z(k,:));
    end
end