findComplexRationalApproximation.m

Go to the documentation of this file.

function [poles,residues,remainder,proportional] = findComplexRationalApproximation(freqBins,freqData)
% This function takes in frequency data and fits a rational approximation
% using Vector Fitting
% freqBins is in Hz
 
% Number of poles to use in fit
numPoles = 12;
 
% Scale incoming frequencies
maxFreq = max(freqBins);
freqBins = freqBins/maxFreq;
 
% Distribute poles evenly across S-plane
poleFreqs = linspace(freqBins(1),freqBins(end),numPoles);
wPoleFreqs = 2*pi*poleFreqs;
sigma = -0.1;
sPoleFreqs = sigma + 1j*wPoleFreqs; %This is a row vector
 
% Or just use real poles
% sPoleFreqs = linspace(-1e9,-0.1e9,numPoles);
 
% Create vector of simulation frequencies
wFreqBins = 2*pi*freqBins;
 
    % Turn this into a matrix, with identical frequencies along the rows
    s = repmat(1j*wFreqBins(:),1,numPoles);
 
for ii = 1:250
    % Ensure poles are stable
    sPoleFreqs = flipReal(sPoleFreqs);
    
    [r,ctilde,Asim] = getResidues(sPoleFreqs,freqData,wFreqBins);
 
    % Extract zeros
    % Form matrix with poles on diagonal
    Az = diag(sPoleFreqs);
    
    Bz = ones(size(sPoleFreqs,2),1);
    
    % Determine zeros (which become new poles)
    z(:,ii) = eig(Az-Bz*ctilde.');
    
    % Set next iteration poles equal to zeros
    sPoleFreqs = z(:,ii).';
end
 
% Ensure poles are stable
sPoleFreqs = flipReal(sPoleFreqs);
 
% Get residues for these final (stable) poles
[r,~,Asim, residual] = getResidues(sPoleFreqs,freqData,wFreqBins);
 
% Extract resiudes: Asim*res = f(s)
freqDataFull = freqData(:);    
r = Asim\freqDataFull;
 
% Compute residual error and print to screen
residual = norm(Asim*r - freqDataFull);
 
% Count number of useful output variables
numVar = numPoles + 2;
 
% Extract these output variables
res = r(1:numVar);
 
%Plot modelled vs actual data
%plot(freqBins,abs(freqData),'-x'); hold on;
%plot(freqBins,abs(Asim*res(:)),'o'); hold off;
 
% Undo scaling of outgoing residues, poles, remainder term and proportional
% terms
poles = sPoleFreqs*maxFreq;
residues = res(1:numVar-2,:)*maxFreq;
remainder = res(numVar-1,:);
proportional = res(numVar,:)/maxFreq;
end