Using tfest to find a system model
Here is some example code that demonstrates how to use 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]
% specify the form of the model to fit. they are in decreasing powers
decreasing powers of s -- numerator for HPF has no constant term, so set it to zero (otherwise tfest will fit a constant value)
%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 -- numerator for HPF has no constant term, so set it to zero (otherwise tfest will fit a constant value) denominatorForm = [ 1 0 ]; % decreasing powers of s -- numerator for HPF has no constant term, so set it to zero (otherwise tfest will fit a constant value) estimatedTransferFunction = FitTransferFunction( magnitude, phase, frequency, numeratorForm, denominatorForm ) 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