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:

%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