Difference between revisions of "Using tfest to find a system model"

The function FitTransferFunction below uses tfest to fit a linear model to frequency response data.

The arguments numeratorForm and denominatorForm are vectors of ones and zeros. Each element corresponds to one term of the numerator or denominator of the transfer function, in order of decreasing powers of $ s $. For example, if you want the function to fit a transfer function of the form $ \frac{K_1 s}{K_2 s + K_3} $ (a high-pass filter), you would use numeratorForm = [ 1 0 ] and denominatorForm = [ 1 1 ]. Using these parameters, prevent tfest will fit a first-order polynomial for the numerator and denominator with no constant term in the numerator. If you wanted to fit a transfer function of the form $ \frac{K_1 s^2}{K_2 s^3 + K_3 s^2 + K_4 s + K_5} $, you would use numeratorForm = [ 1 0 0 ] and denominatorForm = [ 1 1 1 1]

function EstimatedTransferFunction = FitTransferFunction( ...
          MagnitudeRatio, PhaseDifferenceDegrees, FrequencyHertz, NumeratorForm, DenominatorForm )
    % convert magnitude and phase to single complex vector
    complexResponseData = MagnitudeRatio .* exp( 1i .* PhaseDifferenceDegrees .* pi ./ 180 );

    % create optional initial model for tfest so that known zero
    % coefficients can be constrained
    initalModel = idtf( NumeratorForm, DenominatorForm );
    zeroCoefficients = find( NumeratorForm == 0 );
    for ii = 1:numel( zeroCoefficients )
        initalModel.Structure.Numerator.Free(ii) = false;
    end
    zeroCoefficients = find( DenominatorForm == 0 );
    for ii = 1:numel( zeroCoefficients )
        initalModel.Structure.Denominator.Free(ii) = false;
    end

    frequencyResponseData = frd( complexResponseData, FrequencyHertz, 'FrequencyUnit', 'Hz');
    EstimatedTransferFunction = tfest( frequencyResponseData, initalModel );
end

Here is some sample code you can use to test and understand FitTransferFunction.

%Generate some synthetic frequency response measurement data
s = tf( [ 1 0 ], 1 );
highPass = s / ( s + 1 )

[ magnitude, phase, frequency ] = bode( highPass );

% add some random noise
noiseStandardDeviation = 0.05;
magnitude = magnitude + noiseStandardDeviation * randn( size( magnitude ) );
phase = phase + noiseStandardDeviation * randn( size( phase ) );

numeratorForm = [ 1 0 ]; % decreasing powers of s -- first order with no constant term
denominatorForm = [ 1 0 ]; % decreasing powers of s -- first order with constant term

estimatedTransferFunction = FitTransferFunction( magnitude, phase, frequency, numeratorForm, denominatorForm )