getResidues.m

Go to the documentation of this file.

function [residues, ctilde, Asim, residual] = getResidues(poles, freqData, wSimFreqs)
%GETRESIDUES Determines the residues of a rational model fit to frequency
%domain matrix (note that this is true VectorFitting)
%   This function is for use with vector fitting. It simulates a rational
%   model (in adjusted partial fraction form) and determines the resdiues
%   required to fit a function signma(s) using least squares. This can then 
%   in turn be used to update the poles of the rational model for the next 
%   iteration.
    
    % Total number of frequency points
    numFreqs = length(wSimFreqs);
    
    % Number of poles to use in fit
    numPoles = length(poles);
    
    %Turn this into a matrix, with identical frequencies along the rows
    Sk = 1j*repmat(wSimFreqs(:),1,numPoles);
 
    %anorm are the complex starting/current poles
    %Here we create a matrix with identical poles down the columns  
    P = repmat(poles,numFreqs,1);
    
    %Using these poles make 1/(s-a) terms
    %Each row corresponds to a partial fraction model at one frequency
    Asys = 1./(Sk-P);
    
    %To ensure residues come in perfect conjugate pairs, we form A as
    %A = [ANorm + AConj, 1j*ANorm - 1j*AConj];
    
    %How many vectors are we fitting?
    numVec = size(freqData,1)*size(freqData,2);
    numVar = numPoles + 2;
    
    %The partial fraction expansions consits of terms:
    %[1/(s-a1), 1/(s-a2),..., 1, s]
    %This is used as part of the vector fit matrix, but also for evaluating the
    %fit later  
    Asim = [Asys, ones(numFreqs,1), 1j*wSimFreqs(:)*0];
    
    %AsimFull = Asim;
    %freqDataFull = freqData;
    
    %Add the  partial fraction terms for the sigma function 
    %A = [AsimFull, -repmat(Asys,4,1).*repmat(freqDataFull,1,size(Asys,2))];
    A = [Asim, -Asys.*repmat(freqData(:),1,size(Asys,2))];
     
    b = freqData(:);
    
    %Obtain iteration resiudes for model
    x = A\b;
    res = x(1:numVar);
    %res(:,2) = x(1*numVar + 1:2*numVar);
    %res(:,3) = x(2*numVar + 1:3*numVar);
    %res(:,4) = x(3*numVar + 1:4*numVar);
    residues = res(:);
    
    %Determine residues for sigma(s) function
    %ctilde = x(4*numVar + 1:end);
    ctilde = x(end/2+2:end);
    
    residual = norm(A*x-b);
    
    %% Check fit
    %H = simPoleResidueRemainder(wSimFreqs/2*pi, poles, residues);
 
end