Effacer les filtres
Effacer les filtres

How can i multiply the transfer function with a vector)

31 vues (au cours des 30 derniers jours)
Mohamed Nassef
Mohamed Nassef le 20 Nov 2011
How can i multiply the transfer function with a vector ??

Réponse acceptée

Teja Muppirala
Teja Muppirala le 21 Nov 2011
You can use the LSIM command.
For example:
s = tf('s');
G = (1.659e010*s^2 + 9.123e012*s + 4.147e014)/ (1.011e-005*s^7 + 0.02896*s^6 + 99.2*s^5 + 1.461e005*s^4+ 2.187e008*s^3 + 1.042e011*s^2 + 1.308e013*s+ 4.147e014)
t = 0:0.001:10;
u = sin(t.^3);
y = lsim(G,u,t);
plot(t,u,t,y);
legend({'Input (u)', 'Output (y)'})
  1 commentaire
Mohamed Nassef
Mohamed Nassef le 21 Nov 2011
Thanks Teja so much for answering my question. My question is now how can i get the input vector (u) if I have the output vector (y) and the transfer function (G). i mean how can i divide the output vector by the transfer function to get the input again?

Connectez-vous pour commenter.

Plus de réponses (2)

Teja Muppirala
Teja Muppirala le 21 Nov 2011
You need to invert the transfer function. But simply taking 1/G leaves you with more zeros than poles, and so you cannot exactly solve the inverse problem. But you can get a good estimate if you multiply 1/G with a low pass filter to make it proper (number of poles = number of zeros). Then you can call LSIM once again. You will need a high sampling rate to do this accurately.
%%The forward problem
s = tf('s');
G = (1.659e010*s^2 + 9.123e012*s + 4.147e014)/ (1.011e-005*s^7 + 0.02896*s^6 + 99.2*s^5 + 1.461e005*s^4+ 2.187e008*s^3 + 1.042e011*s^2 + 1.308e013*s+ 4.147e014);
t = 0:0.00001:10;
u = sin(t.^3);
y = lsim(G,u,t);
plot(t,u,t,y);
legend({'Input (u)', 'Output (y)'})
%%The inverse problem:
% Step 1. Invert G
G_inverse = 1/G;
% Step 2. Make the low pass filter
lpf_order = numel(pole(G)) - numel(zero(G));
lpf_freq = 10*max(abs(real([pole(G); zero(G)]))); %Cutoff frequency 10x higher than highest zero/pole in G
[num,den] = butter(lpf_order,lpf_freq,'s');
% Step 3. Multiply the inverse by the low pass filter
G_inverse = G_inverse * tf(num,den)
% Step 4. Use LSIM with 'y' as the input
u_estimate = lsim(G_inverse,y,t);
figure;
plot(t,u,t,y,t,u_estimate);
legend({'Actual Input (u)', 'Output (y)', 'Estimated input (u estimate)'})
  1 commentaire
Mohamed Nassef
Mohamed Nassef le 21 Nov 2011
thanks so much..really appreciate it ..it works for me

Connectez-vous pour commenter.


Alex
Alex le 20 Nov 2011
What form is your transfer function in?
  1 commentaire
Mohamed Nassef
Mohamed Nassef le 21 Nov 2011
Dear Alex thanks for your help but it didnt work...the transfer function is in s domain and as follows : TF= (1.659e010 s^2 + 9.123e012 s + 4.147e014)/ (1.011e-005 s^7 + 0.02896 s^6 + 99.2 s^5 + 1.461e005 s^4+ 2.187e008 s^3 + 1.042e011 s^2 + 1.308e013 s+ 4.147e014)
when i apply ilaplace it gives me the following result and inside it there is term called r3 and i dont know what is r3 and still dont know how to multiply the result with a vector in time domain
sum((382493238368367552757760000000000000000*exp(r3*t) + 8414482309222611969638400000000000000*r3*exp(r3*t) + 15301574209142073065472000000000000*r3^2*exp(r3*t))/(65273803904821246875*r3^6 + 160265312512388584439808*r3^5 + 457479253027996880076800000*r3^4 + 539013861833793098219520000000*r3^3 + 605145439338041840762880000000000*r3^2 + 192215073248053527838720000000000000*r3 + 12064170624206046756864000000000000000), r3 in RootOf(s3^7 + (26710885418731430739968*s3^6)/9324829129260178125 + (146393360968959001624576*s3^5)/14919726606816285 + (14373702982234482619187200*s3^4)/994648440454419 + (7172094095858273668300800000*s3^3)/331549480151473 + (30754411719688564454195200000000*s3^2)/2983945321363257 + (1286844866581978320732160000000000*s3)/994648440454419 + 122397836277877616882483200000000000/2983945321363257, s3))

Connectez-vous pour commenter.

Catégories

En savoir plus sur Numerical Integration and Differential Equations dans Help Center et File Exchange

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by