How can ı get transfer function.
65 views (last 30 days)
Show older comments
bahadir safak
on 27 Feb 2018
Commented: Star Strider
on 6 Mar 2018
Dear sir.
I have input data on time and I have output data on time.
Have can I find transfer function. I attached data.csv and graphic.
Best regard.
0 Comments
Accepted Answer
Star Strider
on 27 Feb 2018
This is an interesting problem, and eminently solvable with the System Identification Toolbox.
First, the Fourier transform indicates that there are 3 poles, one at the origin, one at about 550 KHz, and one at infinity (actually the Nyquist frequency, that is for all practical purposes infinity), and 2 zeros, between each pair of the poles. That is all we need to estimate the transfer function.
The Code —
[D,S,R] = xlsread('data.csv');
tv = D(:,1);
vi = D(:,2);
vi(1) = 0;
vo = D(:,3);
[ti,ia,ic] = uniquetol(tv, 4.8E-7); % Get Unique Times
vi = vi(ia); % Corresponding Input Voltages
vo = vo(ia); % Corresponding Output Voltages
Ts = mean(diff(ti)); % Sampling Interval (sec)
ti = [(0 : 5E-7 : 1E-6)'; ti+1.5E-6]; % Pad All Vectors With 3 Initial Time Samples
vi = [zeros(3,1); vi]; % Pad All Vectors With 3 Initial Time Samples
vo = [zeros(3,1); vo]; % Pad All Vectors With 3 Initial Time Samples
Fs = 1/Ts; % Sampling Frequency (Hz)
Fn = Fs/2; % Nyquist Frequency
L = numel(ti); % Vector LEngth
FTvi = fft(vi)/L; % Input Fourier Transform
FTvo = fft(vo)/L; % Output Fourier Transform
FTtf = FTvo./FTvi; % Transfer Function Fourier Transform
Fv = linspace(0,1,fix(L/2)+1)*Fn; % Frequency Vector
Iv = 1:numel(Fv); % Index Vector
figure(1)
plot(Fv, 20*log10(abs(FTtf(Iv))*2))
xlabel('Frequency (Hz)')
ylabel('Amplitude (dB)')
grid
dataobj = iddata(vo,vi,Ts); % Prepare Data For Identification
tfobj = tfest(dataobj, 3, 2); % Create Transfer Function Object
NumTF = tfobj.Numerator; % Transfer Function Numerator
DenTF = tfobj.Denominator; % Transfer Function Numerator
TF = tf(tfobj) % Transfer Function
figure(2)
plot(ti, vi, ti, vo) % Plot Input & Output
xlabel('Time (sec)')
ylabel('Amplitude (V)')
grid
[mag,phase,wout] = bode(tfobj);
figure(3)
subplot(2,1,1)
semilogx(wout, 20*log10(squeeze(mag)), '-b', 'LineWidth',1) % Bode Plot: Magnitude
ylabel('H(f) (dB)')
grid
subplot(2,1,2)
semilogx(wout, squeeze(phase), '-b', 'LineWidth',1) % Bode Plot: Phase
xlabel('Frequency (rad/sec)')
ylabel('Phase (°)')
grid
opt = stepDataOptions;
opt.StepAmplitude = vi(end);
[y,t] = step(tfobj, ti(end), opt); % Calculate Step Response
figure(4)
plot(t,squeeze(y)) % Step Response
title('Step Response of Estimated Transfer Function')
xlabel('Time (sec)')
ylabel('Amplitude (V)')
grid
produces:
TF =
From input "u1" to output "y1":
99.51 s^2 + 1.448e06 s + 3.853e09
--------------------------------------
s^3 + 4612 s^2 + 1.002e07 s + 1.076e10
Continuous-time transfer function.
The transfer function polynomial coefficients with full precision are in ‘NumTF’ and ‘DenTF’.
and this plot:
I will let you refine this, since you know your system better than I do.
6 Comments
More Answers (0)
See Also
Categories
Find more on Data Preparation Basics in Help Center and File Exchange
Products
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!