Obtain transfer function from (noisy) measurement data and use this transfer function for control in simulink

I have a setup with an accelerometer and a speaker. The speaker is used excited with a random noise, and I measure the acceleration in g (after normalization of the output voltage). My question now is what would be the best solution to obtain the transfer function of the given real system in a certain frequency range? Ideally, I would like to avoid to have to guess the number of poles and zeros with a function like tfest.

Currently I have the following code for the estimation of the transfer function. With this solution I was able to obtain the transfer function, but I am at a loss how to implement it into Simulink.

load('datset_y_u.mat')

disp("calculating psd of input/output")
nfft = 32 * samples_per_frame_sensor;
window = nfft;
[pxx,fxx] = pwelch(excitation_signal,window, [], nfft, f_s_audio_sensor);
[pyy,fyy] = pwelch(raw_sensor_signal,window, [], nfft, f_s_audio_sensor);

disp("transfer function estimation")
[Txy,f] = tfestimate(pxx,pyy,window/2,[],nfft/2,f_s_audio_sensor);

figure(1);

subplot(3,1,1);
title("input");
plot(fxx, 10*log10(pxx))
xlabel('Frequency (Hz)')
ylabel('PSD (dB/Hz)')
xlim([0 2000]);

subplot(3,1,2);
title("output");
plot(fyy, mag2db(abs(pyy)))
xlim([0 2000]);
xlabel('Frequency (Hz)')
ylabel('PSD (dB/Hz)')

subplot(3,1,3);
title("transfer function");

plot(f,mag2db(abs(Txy)))
xlim([0 2000]);
xlabel('Frequency (Hz)')
ylabel('PSD (dB/Hz)')

Answer 1

[Sam Chak]

Hi @tiessen

Do you want to estimate the frequency response the Transfer Function or a real mathematical function-type of the Transfer Function?

If it is the latter, then you should consider using the tfest() function to obtain the numerator and denominator coefficients of estimated model. If you know behavior of the system, then you can try guessing the number of poles and zeros.

Having the numerator and denominator coefficients, in Simulink, you can input the parameters in the Transfer Function block.

It is also possible to convert Frequency-Response data into Transfer Function.

% estimate a transfer function Gp with the number of poles and the number of zeros
Gp = tfest(data, np, nz);

% If np < 4, then try using pidtune to autotune the PID gains for the estimated Gp
Gc = pidtune(Gp, 'pid');

[tiessen]

Thank you for your reply :) The final goal is to have a mathmatical representation of the real system (TF) to design a controller. The controller should control the amplitude of the speaker in order to have the same acceleration accross a certain freqeuncy band.

I tried to get the frequency response out of the measurement data and the code in my initial post. For some reason also the transfer function I get with tfestimate and G = Y/U nearly correspond to Y. I wanted to use this data to get a bode plot to approximate the transfer function from that.

Anther question I got is if there is a way to implement the vector I get from tfestimate as a transfer function in simulink?

Answer 2

[Rajiv Singh]

Automatic order determination is not easy/trivial in general. The closest we have to this is automatic selection of state-space model order in ssest/n4sid commands, as in model = ssest(data, 'best')

[tiessen]

Thank you I will look into that.