www.gusucode.com > signal 工具箱matlab源码程序 > signal/@dfilt/@cascade/tosysobj.m

    function Hs = tosysobj(this,returnSysObj)
%TOSYSOBJ Convert to a System object

%   Copyright 2013-2014 The MathWorks, Inc.

    if ~returnSysObj
        % If returnSysObj is false, then it means that we want to know if the
        % System object conversion is supported for the class at hand.

        Hs = true;
        return;
    end

    t =  getfmethod(this);
    coupledAllPassDesign = false;
    if(isprop(t, 'FilterStructure'))
        desString = t.FilterStructure;
        switch(desString)
          case 'cascadeallpass',
            Hs = dsp.CoupledAllpassFilter('Minimum multiplier');
            coupledAllPassDesign = true;
          case 'cascadewdfallpass',
            Hs = dsp.CoupledAllpassFilter('Wave Digital Filter');
            coupledAllPassDesign = true;
        end
    end

    if coupledAllPassDesign
        % Sets properties of dsp.CoupledAllpassFilter matching filter structure in
        % source dfilt.cascade, including
        % - Coefficients (for branches with dfilt.cascade(wdf)allpass
        % - Presence of delays in branches, with various implementations in source
        %   dfilt.cascade, including through dfilt.delay and dfilt.dffir
        % - Individual branch gains, with various implementations in source
        %   dfilt.cascade, including through dfilt.scalar and dfilt.dffir
        setProperties(this, Hs);

        % Generates a warning if source dfilt has persistent memory and non-trivial
        % internal states, since dsp.CoupledAllpassFilter does not support nonzero
        % initial conditions
        warnIfStatesNontrivial(this);
    else
        % Create a cascade of System objects
        if this.nstages == 0
            error(message('signal:dfilt:cascade:tosysobj:noStages'));
        end
        Hs = dsp.FilterCascade(sysobj(this.Stage(1)));
        for idx = 2:this.nstages
            Hs.addStage(sysobj(this.Stage(idx)));
        end
    end


end

% --- Main helper functions ---
function setProperties(dfiltObj, sysObj)

% Parse dfiltObj and get information on branches, delays and gains
    dfiltObj = preprocessCAP(dfiltObj);
    info = getInfoOnAllpassbranches(dfiltObj);

    % Coefficients and delays in Branch #1. Branch1 #1 can be a pure delay
    if(~info(1).IsDelay)
        % Branch not delay - assume info.Branch1.Filter is a
        % dfilt.cascade(wdf)allpass
        % Get coefficients for source dfilt in Branch1 as cell array
        C = getCoefficientsFromSourceDfilt(info(1));
        % If branch also includes delay, add further allpass stage to reflect
        % that. Cdelay is either empty or a cell array
        Cdelay = coefficientValuesForDelay(info(1));
        C = [C, Cdelay];
        % Set constructed Coefficients cell array in target sysObj, Branch #1
        setRelevantCoefficientsInTargetSysobj(C, sysObj, 1);
    else
        % If Branch1 is simple delay, set PureDelayBranch to true and set the
        % Delay property to match the DelayLength found
        sysObj.PureDelayBranch = true;
        sysObj.Delay = info(1).DelayLength;
    end

    % Coefficients and delays for Branch #2 - this must contain a nontrivial
    % cascaded allpass filter
    % Get coefficients for source dfilt in Branch2 as cell array
    C = getCoefficientsFromSourceDfilt(info(2));
    % If branch also includes delay, add further allpass stage to reflect
    % that. Cdelay is either empty or a cell array
    Cdelay = coefficientValuesForDelay(info(2));
    C = [C; Cdelay];

    % Set relevant coefficients in System object (Branch #2), using available
    % coefficients cell array
    setRelevantCoefficientsInTargetSysobj(C, sysObj, 2);

    % Set extra branch gains in target System object
    setBranchGains(info, sysObj);

end
function warnIfStatesNontrivial(dfiltObj)
% This does not do anything, except throwing a warning if source dfilt
% object has populated internal states - Those will not be copied across to
% target System object

% To start, need to get to an actual nontrivial firlter object in the
% structure. Choose to go for the first reachable cascade(wdf)allpass
    branchinfo = getEmptyInfo;

    branchinfo = findFirstAllpassBranch(dfiltObj, branchinfo);

    firstAllpassBranch = extractDfiltAtPath(dfiltObj, branchinfo.Path);

    % Check if warning needed and throw if appropriate
    if firstAllpassBranch.PersistentMemory
        % Get internal states
        DS = getinitialconditions(firstAllpassBranch);
        nontrivialStates = ~isempty(DS) && ~all(all(DS == 0));
        if(nontrivialStates)
            % The conversion to System object ignores
            % non-zero internal states as dsp.CoupledAllpassFilter does not
            % support non-zero initial conditions
            warning(message('signal:dfilt:cascade:tosysobj:nostatesupport'))
        end
    end

    % --- Secondary helper functions ---

end
function emptyinfo = getEmptyInfo

    emptyinfo = struct(...
        'Path', [], ...
        'Filter', [], ...
        'IsDelay', [], ...
        'DelayLength', [], ...
        'Multiplier', []);

end
function info = getInfoOnAllpassbranches(idfilt)
% Parse idfilt and return information on the two parallel
% allpass branches in the structure, including
% - dfilt.cascade(wdf)allpass objects for each branch, where applicable
% - Whether or not each branch is a pure delay block and if yes, the
%   assaciated latency in sampling periods
% - The amound of additional delay to the allpass filter in the branch, if
%   any
% This assumes that
% - The overall structure of idfilt includes a dfilt.parallel
% - That dfilt.parallel has 2 stages
% - At least one of those stages includes (e.g. either directly or within a
%   further dfilt.cascade) a dfilt.cascade(wdf)allpass or a dfilt.dffir
% - If only one of the parallel stages includes a nontrivial dfilt object,
%   then the other includes a dfilt.delay

% Find structural path of up to two cascade(wdf)allpass branches
    info = getPathsOfAllpassBranches(idfilt);

    % If two cascade(wdf)allpass branches are found, return collected info
    if(isempty(info(2).Path))
        % Empty paths mean one allpass not found. Assume
        % - One cascade(wdf)allpass is always found
        % - If only one is found, then a delay is in other branch
        % If only one is found, keep it in Branch2 (make sure only Branch1 can
        % be set to pure delay
        tmp = info(2);
        info(2) = info(1);
        info(1) = tmp;

        % Get delay in Branch1 as parallel with cascade(wdf)allpas in Branch2
        delay1info = getDelayInParallelWithPath(idfilt, info(2));
        info(1).Path = delay1info.Path;
        info(1).IsDelay = true;
    else
        % Get delay in Branch1 a sum of delaying objects in series with the
        % cascade(wdf)allpass filter identified
        delay1info = getDelayInSeriesWithPath(idfilt, info(1));
    end
    info(1).DelayLength = delay1info.DelayLength;

    % Get delay in series with cascade(wdf)allpas in Branch2
    delay2info = getDelayInSeriesWithPath(idfilt, info(2));
    if(delay2info.IsDelay)
        info(2).DelayLength = delay2info.DelayLength;
    end

    % Get gains in series with paths
    gain1info = getGainInSeriesWithPath(idfilt, info(1));
    info(1).Multiplier = gain1info.Multiplier;

    gain2info = getGainInSeriesWithPath(idfilt, info(2));
    info(2).Multiplier = gain2info.Multiplier;

end
function info = getPathsOfAllpassBranches(idfilt)
% Returns [1x2] cell array with paths to up to two non-trivial (of type
% dsp.cascadeallpass, dsp.cascadewdfallpass or dfilt.dffir) leaf filters in
% overall filter structure

    info = repmat(getEmptyInfo, 1, 2);

    info(1) = findFirstAllpassBranch(idfilt, info(1));

    tmp = copy(idfilt);

    tmp = removeCascadeallpassFromPath(tmp, info(1).Path);

    info(2) = findFirstAllpassBranch(tmp, info(2));

end
function oinfo = findFirstAllpassBranch(idfilt, iinfo)
% Returns a structural path pointing to the first cascade(wdf)allpass found
% in the structure hierarchy, as a vector of nested Stage numbers
    oinfo = iinfo;
    ipath = iinfo.Path;
    switch(class(idfilt))
        case {'dfilt.cascade', 'dfilt.parallel'}
            %Try each stage. Stop when we find the first.
            oinfo.Filter = [];
            for ii=1:numel(idfilt.Stage)
                oinfo.Path = [iinfo.Path, ii];
                oinfo = findFirstAllpassBranch(idfilt.Stage(ii), oinfo);
                if ~isempty(oinfo.Path)
                    break; %we found one.
                end
            end    
      case {'dfilt.cascadeallpass', 'dfilt.cascadewdfallpass', 'dfilt.allpass', 'dfilt.wdfallpass'}
        oinfo.Path = ipath;
        oinfo.Filter = idfilt;
        oinfo.IsDelay = false;
      otherwise % including dfilt.scalar, dfilt.dffir, dfilt.delay
        oinfo.Path = [];
        oinfo.Filter = [];

    end

end
function odfilt = removeCascadeallpassFromPath(idfilt, path)
% Places a dfilt.scalar(0) at the location pointed to by the path vector

    if(isempty(path))
        odfilt = idfilt;
        return
    end

    f = idfilt;
    for k = 1:length(path)-1
        f = f.Stage(path(k));
    end
    f.Stage(path(end)) = dfilt.scalar(0);
    odfilt = idfilt;

end
function sergaininfo = getGainInSeriesWithPath(idfilt, mainbranchinfo)

    sergaininfo = getEmptyInfo;

    path = mainbranchinfo.Path;

    branchDfilt = extractParallelBranchContainingPath(idfilt,path);

    gain = getTotalGainInCascade(branchDfilt);
    sergaininfo.Multiplier = gain;

end
function serdelayinfo = getDelayInSeriesWithPath(idfilt, mainbranchinfo)

    serdelayinfo = getEmptyInfo;

    path = mainbranchinfo.Path;

    branchDfilt = extractParallelBranchContainingPath(idfilt, path);

    delaytaps = getTotalDelayInCascade(branchDfilt);

    serdelayinfo.Path = path;
    serdelayinfo.IsDelay = true;
    serdelayinfo.DelayLength = delaytaps;

end
function branchcandidate = extractParallelBranchContainingPath(...
    dfiltObj,path)

    branchcandidate = extractDfiltAtPath(dfiltObj, path);
    for kh = length(path)-1:-1:1
        tmpdfilt = extractDfiltAtPath(dfiltObj, path(1:kh));
        if(isa(tmpdfilt, 'dfilt.parallel'))
            % We hit the top of where we want to search
            break
        end
        branchcandidate = tmpdfilt;
    end

end
function pardelayinfo = getDelayInParallelWithPath(idfilt, mainbranchinfo)

    pardelayinfo = getEmptyInfo;

    path = mainbranchinfo.Path;

    % Find parallel branch by walking upwards from given hierarchy path
    for kh = length(path)-1:-1:1
        pdfilt = extractDfiltAtPath(idfilt, path(1:kh));
        if(isa(pdfilt, 'dfilt.parallel'))
            break
        end
    end
    if(isa(pdfilt, 'dfilt.parallel'))
        curBranch = path(kh + 1);
        parBranch = mod(curBranch,2)+1;
        % Here is the top-level parallel branch to the one passed as input
        parpath = [path(1:kh), parBranch];

        parbranch_dfilt = extractDfiltAtPath(idfilt, parpath);

        delaytaps = getTotalDelayInCascade(parbranch_dfilt);

        pardelayinfo.Path = parpath;
        % pardelayinfo.Filter = parbranch_dfilt;
        pardelayinfo.IsDelay = true;
        pardelayinfo.DelayLength = delaytaps;
    end

end
function odfilt = extractDfiltAtPath(idfilt, path)
    if(~isempty(path))
        odfilt = idfilt.Stage(path(1));
        for k = 2:length(path)
            odfilt = odfilt.Stage(path(k));
        end
    else
        %     odfilt = [];
        odfilt = idfilt;
    end

end
function delaytaps = getTotalDelayInCascade(varargin)
% getTotalDelayInCascade(indfilt, indelay)
    assert(nargin == 1 || nargin == 2)
    indfilt = varargin{1};

    if(nargin == 2)
        indelay = varargin{2};
    else
        indelay = 0;
    end

    switch(class(indfilt))
      case 'dfilt.cascade',
        % Sum delay over all cascade stages
        numstages = nstages(indfilt);
        delaytaps = indelay;
        for k = 1:numstages
            delaytaps = getTotalDelayInCascade(indfilt.Stage(k), delaytaps);
        end
      case 'dfilt.parallel',
        % Assuming to only be dealing with a cascade of dfilt.stages
        assert(false)
      case 'dfilt.scalar',
        delaytaps = indelay;
      case 'dfilt.delay',
        delaytaps = indelay + indfilt.Latency;
      case 'dfilt.dffir',
        delaytaps = indelay + median(grpdelay(indfilt, 3));
      otherwise
        % Assume no delay to report for other filter classes
        delaytaps = indelay;
    end

end
function gain = getTotalGainInCascade(varargin)
% getTotalGainInCascade(indfilt, ingain)
    assert(nargin == 1 || nargin == 2)
    indfilt = varargin{1};

    if(nargin == 2)
        ingain = varargin{2};
    else
        ingain = 1;
    end

    switch(class(indfilt))
      case 'dfilt.cascade',
        % Sum delay over all cascade stages
        numstages = nstages(indfilt);
        gain = ingain;
        for k = 1:numstages
            gain = getTotalGainInCascade(indfilt.Stage(k), gain);
        end
      case 'dfilt.parallel',
        % Assuming to only be dealing with a cascade of dfilt.stages
        assert(false)
      case 'dfilt.scalar',
        gain = ingain * indfilt.Gain;
      case 'dfilt.delay',
        gain = ingain;
      case 'dfilt.dffir',
        % Only report gain if this has Numerator of type [0 ... 0 gain]
        if(all(indfilt.Numerator(1:end-1) == 0))
            gain = ingain * indfilt.Numerator(end);
        else
            gain = ingain;
        end
      otherwise
        % Assume no additional gain to report
        gain = ingain;
    end

end
function C = getCoefficientsFromSourceDfilt(sourceBranchInfo)

    CascadeFilter = sourceBranchInfo.Filter;
    % Number of filter sections in top branch
    c = CascadeFilter.AllpassCoefficients;
    if isstruct(c)
        numsec = length(fieldnames(c));
    else
        numsec = 1; %dfilt.allpass case
    end
    C = cell(numsec, 1);
    for k = 1:numsec
        if isa(CascadeFilter, 'dfilt.allpass') || isa(CascadeFilter, 'dfilt.wdfallpass')
            C{k} = CascadeFilter.AllpassCoefficients;
        else    
            C{k} = CascadeFilter.AllpassCoefficients.(['Section',num2str(k)]);
        end
        if(isa(CascadeFilter, 'dfilt.cascadewdfallpass'))
            C{k} = dsp.AllpassFilter.poly2wdf(C{k});
        end
    end

end
function setRelevantCoefficientsInTargetSysobj(...
    Coefficients, sysObj, targetBranchNumber)

    switch(sysObj.Structure)
      case 'Minimum multiplier'
        sysObj.(['AllpassCoefficients', num2str(targetBranchNumber)]) = ...
            Coefficients;
      case 'Wave Digital Filter'
        sysObj.(['WDFCoefficients', num2str(targetBranchNumber)]) = ...
            Coefficients;
      otherwise
        % Only filter classes dfilt.cascadeallpass and
        % dfilt.cascadewdfallpass expected
        assert(false)
    end

end
function DelayCoefficients = ...
        coefficientValuesForDelay(sourceBranchInfo)

    if(sourceBranchInfo.DelayLength > 0)
        if(isa(sourceBranchInfo.Filter, 'dfilt.cascadeallpass'))
            % The task is easy in this case as the 'Minimum multiplier'
            % Structure in dsp.AllpassFilter and dsp.CoupledAllpassFilter does
            % support arbitrary long single sections. Simply use one section
            DelayCoefficients = {zeros(1, sourceBranchInfo.DelayLength)};

        elseif(isa(sourceBranchInfo.Filter, 'dfilt.cascadewdfallpass'))
            % Things are more complicated in this case - the 'Wave Digital
            % Filter' Structure in dsp.AllpassFilter and
            % dsp.CoupledAllpassFilter only supports sections of order 1, 2 and
            % 4. Use as many 2nd-order sections as necessary plus possibly one
            % first-order section. All inner coefficients to be zeros

            % First take care of teh 2nd-order sections
            nSecondOrderSections = floor(sourceBranchInfo.DelayLength/2);
            DelayCoefficients = repmat({[0, 0]}, nSecondOrderSections, 1);
            % Add last 1st order section if necessary
            if(sourceBranchInfo.DelayLength - 2*nSecondOrderSections == 1)
                DelayCoefficients = [DelayCoefficients; {0}];
            end
        else
            % Only calsses dfilt.cascadeallpas and dfilt.cascadewdfallpas
            % expected as inputs
            assert(false)
        end
    else
        DelayCoefficients = [];
    end

end
function setBranchGains(info, sysObj)

    switch(info(1).Multiplier)
      case 1,
        sysObj.Gain1 = '1';
      case -1,
        sysObj.Gain1 = '-1';
      case 1i,
        sysObj.Gain1 = '0+1i';
      case -1i,
        sysObj.Gain1 = '0-1i';
      otherwise
        error(message('signal:dfilt:cascade:tosysobj:badbranchgain'));
    end
    switch(info(2).Multiplier)
      case 1,
        sysObj.Gain2 = '1';
      case -1,
        sysObj.Gain2 = '-1';
      case 1i,
        sysObj.Gain2 = '0+1i';
      case -1i,
        sysObj.Gain2 = '0-1i';
      otherwise
        error(message('signal:dfilt:cascade:tosysobj:badbranchgain'));
    end
end