How can ı get transfer function.

65 views (last 30 days)
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.

Accepted Answer

Star Strider
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
bahadir safak
bahadir safak on 6 Mar 2018
Edited: bahadir safak on 6 Mar 2018
İt works ; Thank you
Star Strider
Star Strider on 6 Mar 2018
As always, my pleasure.

Sign in to comment.

More Answers (0)

Community Treasure Hunt

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

Start Hunting!