Transfer function of the highly distorted system; impact hammer testing and frequency response of the system

Hello, everyone!

For couple of weeks now, I am dealing with transfer function estimation. The system which I am currently examining is highly distorted at high frequencies due to its structural dynamics. I analysed its frequency response and it is getting really messy above 2,5 kHz.

I want to compensate this system with the use of the transfer function. I have some input and output data (attached you will find the graph illustrating input and output signals). It has been carried out by means of the impact hammer testing. To cut the long story short: input - hammer hit, output - response of the system (oscillations).

I tried to handle it with the System Identification Toolbox. What I got so far is:

Could you please provide me with any hints how to get this function into time domain? The aim of my work is to reproduce the input signal without the distortions.

Thank you very much for any feedback.

Answer 1

[Raj]

Not sure why you took the trouble of posting so much details. From what I understand, the one line question is that you have got a transfer function in 'S' domain and you want to convert this to time domain. You can use the ilaplace command to do that.

Note: You should have symbolic math toolbox in order to use this.

[Daniel]

Thank you for your answer, Raj!

I was trying to do it with "ilaplace" command. However, what I get then is:

val =
 
-(16195*symsum((exp(t*root(s3^3 + (12249*s3^2)/25000 + (19438*s3)/5 + 177760000, s3, k))*root(s3^3 + (12249*s3^2)/25000 + (19438*s3)/5 + 177760000, s3, k)^2)/(2*(37500*root(s3^3 + (12249*s3^2)/25000 + (19438*s3)/5 + 177760000, s3, k)^2 + 12249*root(s3^3 + (12249*s3^2)/25000 + (19438*s3)/5 + 177760000, s3, k) + 48595000)), k, 1, 3))/2 - 6076750000000*symsum(exp(t*root(s3^3 + (12249*s3^2)/25000 + (19438*s3)/5 + 177760000, s3, k))/(2*(37500*root(s3^3 + (12249*s3^2)/25000 + (19438*s3)/5 + 177760000, s3, k)^2 + 12249*root(s3^3 + (12249*s3^2)/25000 + (19438*s3)/5 + 177760000, s3, k) + 48595000)), k, 1, 3) - 70327500*symsum((exp(root(s3^3 + (12249*s3^2)/25000 + (19438*s3)/5 + 177760000, s3, k)*t)*root(s3^3 + (12249*s3^2)/25000 + (19438*s3)/5 + 177760000, s3, k))/(2*(12249*root(s3^3 + (12249*s3^2)/25000 + (19438*s3)/5 + 177760000, s3, k) + 37500*root(s3^3 + (12249*s3^2)/25000 + (19438*s3)/5 + 177760000, s3, k)^2 + 48595000)), k, 1, 3)

And this is the part I am struggling with. Do you happen to know the reason?

[Raj]

Can you post your code? Or atleast the part where you are taking inverse Laplace of the 'S' domain transfer function?

[Daniel]

Of course!

% Transfer function estimation
data = iddata(output_data, input_data, 0.00001); % Sampling rate was 100 kHz
sys = tfest(data, 3);

% Inverse Laplace transform
syms s
tf_s =  (-3.24*10^4*s^2 - 2.811*10^8*s - 2.423*10^13) / (s^3 + 4.879*10^4*s^2 + 3.866*10^8*s + 1.77*10^13)
tf_time = ilaplace(tf_s)

As a results, I get:

val =
 
- 281100000*symsum((exp(root(s3^3 + 48790*s3^2 + 386600000*s3 + 17700000000000, s3, k)*t)*root(s3^3 + 48790*s3^2 + 386600000*s3 + 17700000000000, s3, k))/(97580*root(s3^3 + 48790*s3^2 + 386600000*s3 + 17700000000000, s3, k) + 3*root(s3^3 + 48790*s3^2 + 386600000*s3 + 17700000000000, s3, k)^2 + 386600000), k, 1, 3) - 32400*symsum((exp(t*root(s3^3 + 48790*s3^2 + 386600000*s3 + 17700000000000, s3, k))*root(s3^3 + 48790*s3^2 + 386600000*s3 + 17700000000000, s3, k)^2)/(3*root(s3^3 + 48790*s3^2 + 386600000*s3 + 17700000000000, s3, k)^2 + 97580*root(s3^3 + 48790*s3^2 + 386600000*s3 + 17700000000000, s3, k) + 386600000), k, 1, 3) - 24230000000000*symsum(exp(t*root(s3^3 + 48790*s3^2 + 386600000*s3 + 17700000000000, s3, k))/(3*root(s3^3 + 48790*s3^2 + 386600000*s3 + 17700000000000, s3, k)^2 + 97580*root(s3^3 + 48790*s3^2 + 386600000*s3 + 17700000000000, s3, k) + 386600000), k, 1, 3)

As you can see it above, I decided to rewrite the function manually as I did not know how to get numerator and denumerator in the correct data type. Apart from that, please have a look at it.

Thank you for your help.

[Star Strider]

If you add:

Val = vpa(val)

you get:

Val =
- 32193.426574018597179244606329801*exp(-48363.599240064750874320653505083*t) - 206.57342598140282075539367019915*exp(-213.20037996762456283967324745842*t)*cos(19129.356239267807979944343311528*t) - 13452.518699763963342500280154929*exp(-213.20037996762456283967324745842*t)*sin(19129.356239267807979944343311528*t)

that in a slightly more readable form is:

[Daniel]

Thank you, it worked! But another question is how to use this transfer function to correct the data in order to reproduce undistorted input signal? As transfer function is represented by the equation:

I just thought that I would get it simply by dividing the output signal by the transfer function in time domain. Is my train of thought correct?

[Raj]

"dividing the output signal by the transfer function in time domain" - No that is not correct! Transfer function is output/input in 'S' domain not in time domain. In time domain 'convolution' will come into picture. See this for better understanding:

[Daniel]

close all
clear all
clc

t = linspace(0,2,200000); % Sampling rate: 100 kHz; t - time [s]
t = t';

% Import the data and decrease the number of samples
input = importdata('C:\Users\Lechowicz\Desktop\input.txt', '');
output = importdata('C:\Users\Lechowicz\Desktop\output.txt', '');

% Plot the data
figure('Name', 'Original Data', 'NumberTitle', 'on');
subplot(2,1,1)
plot(t, input)
title('Input Signal');
xlabel('Time [s]');
ylabel('Force [N]');
subplot(2,1,2)
plot(t, output)
title('Output Signal');
xlabel('Time [s]');
ylabel('Force [N]');

% Estimate the transfer function and convert it to time domain
data = iddata(output, input, 0.00001); 
sys = tfest(data, 3, 2);
[num, den] = tfdata(sys, 'v');
syms s
tf_s = poly2sym((num),s)/poly2sym((den),s);

% Get the transfer function into time domain
time_domain_tf = ilaplace(tf_s);
time_domain_tf = vpa(time_domain_tf);

% Create matrix of time domain transfer function
for i = (linspace(0,2,200000))'
    tf_time_domain_matrix = 10.753107173906815697416094377728.*exp(-3.327484297190415992653272529398.*i) - 840.46923468767926446235542117068.*exp(-246.19890346395301270080891184263.*i).*cos(19126.216965822094847808546864243.*i) - 14175.620532280572917636521940787.*exp(-246.19890346395301270080891184263.*i).*sin(19126.216965822094847808546864243.*i);
end

% I am almost sure I messed up here
output_corrected = conv(tf_time_domain_matrix, output);

% Plot the data
figure('Name', 'Reproduced Input Signal', 'NumberTitle', 'on');
plot(linspace(0,2,length(output_corrected)), output_corrected);
title('Reproduced Output Signal');
xlabel('Time [s]');
ylabel('Force [N]');

I give up. Could you please point out where I did mistake?

As far as I understand it, by correcting my output data I should get something that looks like my input data, i.e. single impact. I can't achieve it for nuts.

I attached my data as well.

All in all, it's getting better and better! Thank you very much for your help.