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.
