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. The model is a transfer function whose numerator and denominator are polynomials in $ s $. (Recall that $ s=j\omega $.) Each element of the vectors corresponds to a single term of the numerator or denominator, 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 ]. These parameters, 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 = zeroCoefficients
        initalModel.Structure.Numerator.Free(ii) = false;
    end
    zeroCoefficients = find( DenominatorForm == 0 );
    for ii = 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, omega ] = bode( highPass );
frequency = omega / ( 2 * pi );

% 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 1 ]; % decreasing powers of s -- first order with constant term

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