classdef (CaseInsensitiveProperties=true, TruncatedProperties=true) sos < matlab.mixin.SetGet & matlab.mixin.Copyable
%lwdfilt.sos class
%    lwdfilt.sos properties:
%       sosMatrix - Property is of type 'mxArray'  
%       ScaleValues - Property is of type 'mxArray'  
%       refsosMatrix - Property is of type 'mxArray'  
%       refScaleValues - Property is of type 'mxArray'  
%
%    lwdfilt.sos methods:
%       computeimpz -   Compute impulse response
%       computephasedelay -   Calculate the phase delay of the filter.
%       computephasez -   Calculate the phase response.
%       impzlength - Length of the impulse response for a digital filter.
%       reffilter -   Returns the double representation of the filter object.
%       thiscoefficients - Filter coefficients.
%       thistf -   Return the transfer function.
%       zpk -  Discrete-time filter zero-pole-gain conversion.


properties (SetObservable)

    %SOSMATRIX Property is of type 'mxArray'
    sosMatrix = [ 1, 0, 0, 1, 0, 0 ];
    %SCALEVALUES Property is of type 'mxArray'
    ScaleValues = [ 1, 1 ]';
    %REFSOSMATRIX Property is of type 'mxArray'
    refsosMatrix = [ 1, 0, 0, 1, 0, 0 ];
    %REFSCALEVALUES Property is of type 'mxArray'
    refScaleValues = [ 1, 1 ]';
end

methods  % constructor block
    function this = sos(sosMatrix, scales)
        %SOS   Construct a SOS object.

        if nargin > 0
            this.sosMatrix = sosMatrix;
            this.refsosMatrix = sosMatrix;
            if nargin > 1
                this.ScaleValues = scales;
                this.refScaleValues = scales;
            end
        end

    end  % sos

    %----------------------------------------
    function [I, T] = computeimpz(this, varargin)
        %COMPUTEIMPZ   Compute impulse response
        % [I, T] = computeimpz(this, N, Fs)

        % Get SOS matrix with embedded scale values
        sosMatrix = getsosmatrix(this); %#ok<*PROPLC>

        [I, T] = impz(sosMatrix,varargin{:});
    end  % computeimpz

    %----------------------------------------
    function [Ph, W] = computephasedelay(this,N,varargin)
        %COMPUTEPHASEDELAY   Calculate the phase delay of the filter.

        if nargin < 2
            N = 8192;
        end

        % Get SOS matrix, set embedScaleValues falg to false
        sosMatrix = getsosmatrix(this,false);

        [Ph, W] = phasedelay(sosMatrix, N, varargin{:});

    end  % computephasedelay

    %----------------------------------------
    function [Ph, w] = computephasez(this, N, varargin)
        %COMPUTEPHASEZ   Calculate the phase response.

        if nargin < 2
            N = 8192;
        end

        % Get SOS matrix, set embedScaleValues falg to false
        sosMatrix = getsosmatrix(this,false);

        [Ph, w] = phasez(sosMatrix, N, varargin{:});

    end  % computephasez

    %----------------------------------------
    function len = impzlength(this, varargin)
        %IMPZLENGTH Length of the impulse response for a digital filter.

        % Get SOS matrix, set embedScaleValues falg to false
        sosMatrix = getsosmatrix(this,false);

        len = impzlength(sosMatrix,varargin{:});

    end  % impzlength

    %----------------------------------------
    function Hd = reffilter(this)
        %REFFILTER   Returns the double representation of the filter object.

        Hd = lwdfilt.sos;
        Hd.sosMatrix = this.refsosMatrix;
        Hd.refsosMatrix = this.refsosMatrix;
        Hd.ScaleValues = this.refScaleValues;
        Hd.refScaleValues = this.refScaleValues;

    end  % reffilter

    %----------------------------------------
    function c = thiscoefficients(Hd)
        %THISCOEFFICIENTS Filter coefficients.
        %   C = THISCOEFFICIENTS(Hd) returns a cell array of coefficients of
        %   discrete-time filter Hd.

        c = {Hd.sosMatrix, Hd.ScaleValues};

    end  % thiscoefficients

    %----------------------------------------
    function [b, a] = thistf(this)
        %THISTF   Return the transfer function.

        sv    = get(this, 'ScaleValues');
        sosm  = get(this, 'sosMatrix');
        nsecs = size(sosm, 1);

        % Pad the scale values with ones
        sv = [sv;ones(nsecs+1-length(sv),1)];

        % Embed gains in matrix
        for n = 1:nsecs
            sosm(n,1:3) = sosm(n,1:3)*sv(n);
        end

        b = 1; a = 1;
        if nsecs>0
            [b, a] = sos2tf(sosm,sv(nsecs+1));
        end

    end  % thistf

    %----------------------------------------
    function [z,p,k] = zpk(this)
        %ZPK  Discrete-time filter zero-pole-gain conversion.
        %   [Z,P,K] = ZP(this) returns the zeros, poles, and gain corresponding to the
        %   discrete-time filter this in vectors Z, P, and scalar K respectively.

        sosm = get(this, 'sosMatrix');

        nsecs = size(sosm, 1);

        if length(this.ScaleValues) > nsecs + 1
            warning(message('signal:lwdfilt:sos:zpk:scalevalues', nsecs + 1));
        end

        z = [];
        p = [];
        k = prod(this.ScaleValues);
        for indx = 1:size(sosm,1)
            [z1,p1,k1] = sos2zp(sosm(indx,:));
            z = [z;z1];
            p = [p;p1];
            k = k*k1;
        end

    end  % zpk

    %----------------------------------------
    function [h, w] = computefreqz(this, N, varargin)
        %COMPUTEFREQZ

        if nargin < 2
            N = 8192;
        end

        % Get SOS matrix with embedded scale values
        sosMatrix = getsosmatrix(this);

        [h, w] = freqz(sosMatrix, N, varargin{:});

    end  % computefreqz

    %----------------------------------------
    function [Gd, w] = computegrpdelay(this, N, varargin)
        %COMPUTEGRPDELAY

        if nargin < 2
            N = 8192;
        end

        % Get SOS matrix, set embedScaleValues falg to false
        sosMatrix = getsosmatrix(this,false);

        [Gd, w] = grpdelay(sosMatrix, N, varargin{:});

    end  % computegrpdelay

    %----------------------------------------
    function [S, T] = computestepz(this, varargin)
        %COMPUTESTEPZ

        % Get SOS matrix with embedded scale values
        sosMatrix = getsosmatrix(this);

        [S, T] = stepz(sosMatrix,varargin{:});

    end  % computestepz

    %----------------------------------------
    function [hz, wz, phiz, opts] = computezerophase(this, N, varargin)
        %COMPUTEZEROPHASE

        % This should be private

        if nargin < 2
            N = 8192;
        end

        % Get SOS matrix with embedded scale values
        sosMatrix = getsosmatrix(this);

        [hz,wz,phiz,opts] = zerophase(sosMatrix, N, varargin{:});

    end  % computezerophase

    %----------------------------------------
    function sosMatrix = getsosmatrix(this,embedScaleValues)
        %GETSOSMATRIX

        if nargin == 1
            embedScaleValues = true;
        end

        sosMatrix = get(this, 'sosMatrix');
        if embedScaleValues
            scales    = get(this, 'ScaleValues');
            sosMatrix = embedsosscalevalues(sosMatrix, scales);
        end
        if size(sosMatrix,1) == 1
            %Analysis only supports SOS matrices with 2 or more sections
            sosMatrix = [sosMatrix ;  [1 0 0 1 0 0]];
        end
    end  % getsosmatrix

end 
end  % classdef

%--------------------------------------------------------------------------
function sos = embedsosscalevalues(sos,g)

if size(sos,2) ~= 6
    error(message('signal:lwdfilt:sos:getsosmatrix:InvalidSOSMatrix'));
end
numSecs = size(sos,1);
if length(g)~= 1 && (length(g) > numSecs+1 || length(g) < numSecs)
    error(message('signal:lwdfilt:sos:getsosmatrix:InvalidScaleValues'));
end

if length(g) == 1
    sos(1,1:3) = sos(1,1:3).*g;
else
    if (length(g) == numSecs+1)
        lastGain = 1;
        if g(end) ~= 1
            % Spread the last non-unity gain across all sections
            lastGain = g(end)^(1/numSecs);
        end
        g = g(1:end-1)*lastGain;
    end
    g = g(:,[1 1 1]);
    sos(:,1:3) = sos(:,1:3).*g;
end
end  % embedsosscalevalues
