Help with tfest function. Resulting transfer function deviates from input on either side of freq range.
Show older comments
Requesting help with tfest options to resolve output transfer function issue where the middle of the range matches input values well, but deviates at both low and high ends of the input frequency range. Could this be due to optional inputs to tfest command?
code:
clear
tic
% create logspaced frequency vector
Fvec = logspace( log10(0.001), log10(1000), 10000 ); % 0.1 - 250 hz
% Equations/Constants to Calculate Ci Gain from frequency vector
% % Lee Pradko Vertical (Z) constants and Gain equations:
K0 = 4.35373;
K1 = 1.356;
Wi = 2*pi*Fvec;
F1 = -0.10245296*10^-9*Wi.^6 + 0.17583343*10^-5*Wi.^4 - 0.44600722*10^-2*Wi.^2 + 1;
F2 = +0.12881887*10^-7*Wi.^4 - 0.93394367*10^-4*Wi.^2 + 0.10543059;
F3 = -0.45416156*10^-9*Wi.^6 + 0.37667129*10^-5*Wi.^4 - 0.56104406*10^-2*Wi.^2 + 1;
F4 = -0.21179193*10^-11*Wi.^6 + 0.5172811*10^-7*Wi.^4 - 0.17946748*10^-3*Wi.^2 + 0.10543059;
Ci = K1.*K0.*(F1.*F4 - F2.*F3) ./ (F3.^2 + Wi.^2.*F4.^2);
sin_phi = (F2.*F3 - F1.*F4).*Wi ./ sqrt( (Wi.^2.*F4.^2+F3.^2).*(F1.^2+Wi.^2.*F2.^2));
phi = asin( sin_phi );
phase = phi* 180/pi;
%create frequency response data
data = frd(Ci.*exp(1j * phi ), Wi);
% Assuming 'data' is your iddata object and 'targetFit' is your desired percentage
targetFit = 99;
np = 1; % Start with 1 pole
nz = 0; % Start with 0 zeros
while true
% Estimate transfer function with current np and nz
sys = tfest(data, np, nz);
% Retrieve the Fit percentage (Normalized RMSE)
currentFit = sys.Report.Fit.FitPercent;
fprintf('Poles: %d, Zeros: %d | Current Fit: %.2f%%\n', np, nz, currentFit);
% Check if target fit is reached
if currentFit >= targetFit
fprintf('Desired fit of %.2f%% reached! Stopping.\n', targetFit);
break;
end
% Update your logic for iterating poles and zeros here.
% Example: Increment zeros up to np-1, then increment poles.
if nz < max(np - 1, 0)
nz = nz + 1;
else
np = np + 1;
nz = 0; % Reset zeros for the new pole count
end
% Set a safe break condition to avoid infinite loops
if np > 15
fprintf('Maximum iterations reached without meeting target fit.\n');
break;
end
end
sys
sys_zpk = zpk(sys)
[mag_raw, phase_raw, wout] = bode(sys,Wi);
mag = squeeze(mag_raw);
phase = squeeze(phase_raw);
Frequency_Hz = wout / (2 * pi);
bodeTable = table(wout, Frequency_Hz, mag, phase);
%writetable(bodeTable, 'BodeDataLeePradkotest.xlsx');
wout_Hz = wout / (2 * pi);
% interpolated magnitude
MAGvec = interp1( wout_Hz, mag, Fvec, 'pchip' );
title_str = 'Magnitude Gain, Z-Vertical (Up-Down) Axis';
lim_y = 15;
% Overlay both curves
set( groot, 'Units', 'inches' )
screensize = get(groot, 'ScreenSize' );
FD = min( screensize(3), screensize(4) );
f1 = figure( 'Units', 'inches', 'Position', [0, 0, FD+1.5, FD] );
movegui( 'center' )
th = get( f1, 'Theme' );
subplot( 2, 1, 1 )
plot( Fvec, Ci, 'r*' )
hold on
plot( Fvec, mag, 'b', 'LineWidth', 3 )
%plot( Fvec, MAGvecL, 'g', 'LineWidth', 3 )
grid on
title( title_str )
xlabel( 'Frequency (Hz)' )
ylabel( 'Gain ( Watts/(ft/s²)² )')
legend( 'Input Data', 'TF Output' )
xlim( [ 0 lim_y ] )
subplot( 2, 1, 2 )
loglog( Fvec, Ci, 'r*' )
hold on
loglog( Fvec, MAGvec, 'b', 'LineWidth', 3 )
%loglog( Fvec, MAGvecL, 'g', 'LineWidth', 3 )
grid on
xlabel( 'Frequency (Hz)' )
ylabel( 'Gain ( Watts/(ft/s²)² )')
legend( 'Input Data', 'TF Output')
toc
Accepted Answer
Hi Jared,
Maybe look into the WeightingFilter opton of tfestOptions. I tried with 'inv' below and seemed to get a result more in line with what's desired, but it doesn't reach the TargetFit (I don't know what that meansn nor if the TargetFit interacts with the WeightingFilter) and there was some other warning, but the results look pretty good. Also, I didn't undersand why interp1 was used, so I got rid of that in the code below, but the interpolated magnitude did look fine when I checked.
clear
tic
% create logspaced frequency vector
Fvec = logspace( log10(0.001), log10(1000), 10000 ); % 0.1 - 250 hz
% Equations/Constants to Calculate Ci Gain from frequency vector
% % Lee Pradko Vertical (Z) constants and Gain equations:
K0 = 4.35373;
K1 = 1.356;
Wi = 2*pi*Fvec;
F1 = -0.10245296*10^-9*Wi.^6 + 0.17583343*10^-5*Wi.^4 - 0.44600722*10^-2*Wi.^2 + 1;
F2 = +0.12881887*10^-7*Wi.^4 - 0.93394367*10^-4*Wi.^2 + 0.10543059;
F3 = -0.45416156*10^-9*Wi.^6 + 0.37667129*10^-5*Wi.^4 - 0.56104406*10^-2*Wi.^2 + 1;
F4 = -0.21179193*10^-11*Wi.^6 + 0.5172811*10^-7*Wi.^4 - 0.17946748*10^-3*Wi.^2 + 0.10543059;
Ci = K1.*K0.*(F1.*F4 - F2.*F3) ./ (F3.^2 + Wi.^2.*F4.^2);
sin_phi = (F2.*F3 - F1.*F4).*Wi ./ sqrt( (Wi.^2.*F4.^2+F3.^2).*(F1.^2+Wi.^2.*F2.^2));
phi = asin( sin_phi );
phase = phi* 180/pi;
%create frequency response data
data = frd(Ci.*exp(1j * phi ), Wi);
% Assuming 'data' is your iddata object and 'targetFit' is your desired percentage
targetFit = 99;
np = 1; % Start with 1 pole
nz = 0; % Start with 0 zeros
while true
% Estimate transfer function with current np and nz
sys = tfest(data, np, nz,tfestOptions('WeightingFilter','inv'));
% Retrieve the Fit percentage (Normalized RMSE)
currentFit = sys.Report.Fit.FitPercent;
%fprintf('Poles: %d, Zeros: %d | Current Fit: %.2f%%\n', np, nz, currentFit);
% Check if target fit is reached
if currentFit >= targetFit
fprintf('Desired fit of %.2f%% reached! Stopping.\n', targetFit);
break;
end
% Update your logic for iterating poles and zeros here.
% Example: Increment zeros up to np-1, then increment poles.
if nz < max(np - 1, 0)
nz = nz + 1;
else
np = np + 1;
nz = 0; % Reset zeros for the new pole count
end
% Set a safe break condition to avoid infinite loops
if np > 15
fprintf('Maximum iterations reached without meeting target fit.\n');
break;
end
end
sys
%sys_zpk = zpk(sys)
[mag_raw, phase_raw, wout] = bode(sys,Wi);
mag = squeeze(mag_raw);
phase = squeeze(phase_raw);
Frequency_Hz = wout / (2 * pi);
wout_Hz = wout / (2 * pi);
% interpolated magnitude
%MAGvec = interp1( wout_Hz, mag, Fvec, 'pchip' );
title_str = 'Magnitude Gain, Z-Vertical (Up-Down) Axis';
lim_y = 15;
% Overlay both curves
figure
subplot( 2, 1, 1 )
plot( Fvec, Ci, 'r*' )
hold on
plot( Fvec, mag, 'b', 'LineWidth', 3 )
%plot( Fvec, MAGvecL, 'g', 'LineWidth', 3 )
grid on
title( title_str )
xlabel( 'Frequency (Hz)' )
ylabel( 'Gain ( Watts/(ft/s²)² )')
legend( 'Input Data', 'TF Output' )
xlim( [ 0 lim_y ] )
subplot( 2, 1, 2 )
loglog( Fvec, Ci, 'r*' )
hold on
%loglog( Fvec, MAGvec, 'b', 'LineWidth', 3 )
loglog( Fvec, mag, 'b', 'LineWidth', 3 )
%loglog( Fvec, MAGvecL, 'g', 'LineWidth', 3 )
grid on
xlabel( 'Frequency (Hz)' )
ylabel( 'Gain ( Watts/(ft/s²)² )')
legend( 'Input Data', 'TF Output')
toc
More Answers (0)
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!
