How to find the Transfer function of a Simulink output plot?

1 view (last 30 days)
Hello My friends
I have a lot of figures and I want to find the transfer functions as I=Imax cos(ωt+ϴ) and V=Vmax cos(ωt+ϴ)
here is the values of V and I
Can you Advise?
I attached one of my plots.
I really appreciate your assistance.

Accepted Answer

Sam Chak
Sam Chak on 28 Sep 2023
Think you want to apply the curve-fitting method to fit the sinusoidal model to the time-series data.
load('V_I.mat')
Vtime = Vabc3.Time;
Vdat1 = Vabc3.Data(:,1);
ta = Vtime(1:1001);
Va = Vdat1(1:1001);
tb = Vtime(1002:2100);
Vb = Vdat1(1002:2100);
tc = Vtime(2101:3001);
Vc = Vdat1(2101:3001);
figure(1)
plot(Vtime, Vdat1), ylim([-1.5 1.5]), grid on
xlabel('Time (seconds)'), ylabel('V_{abc}')
[Vafit, gof1] = fit(ta, Va, 'sin1')
Vafit =
General model Sin1: Vafit(x) = a1*sin(b1*x+c1) Coefficients (with 95% confidence bounds): a1 = 0.998 (0.998, 0.998) b1 = 314.2 (314.2, 314.2) c1 = 1.031 (1.031, 1.031)
gof1 = struct with fields:
sse: 4.5659e-13 rsquare: 1.0000 dfe: 998 adjrsquare: 1.0000 rmse: 2.1389e-08
figure(2)
plot(Vafit, ta, Va), ylim([-1.5 1.5]), grid on
xlabel('Time (seconds)'), ylabel('V_{a}')
[Vbfit, gof1] = fit(tb, Vb, 'sin1')
Vbfit =
General model Sin1: Vbfit(x) = a1*sin(b1*x+c1) Coefficients (with 95% confidence bounds): a1 = 0.1984 (0.1983, 0.1985) b1 = 314.2 (314.1, 314.2) c1 = 7.268 (7.266, 7.27)
gof1 = struct with fields:
sse: 8.6937e-04 rsquare: 1.0000 dfe: 1096 adjrsquare: 1.0000 rmse: 8.9063e-04
figure(3)
plot(Vbfit, tb, Vb), ylim([-0.3 0.3]), grid on
xlabel('Time (seconds)'), ylabel('V_{b}')
[Vcfit, gof1] = fit(tc, Vc, 'sin1')
Vcfit =
General model Sin1: Vcfit(x) = a1*sin(b1*x+c1) Coefficients (with 95% confidence bounds): a1 = 0.9971 (0.9958, 0.9984) b1 = 314.2 (314.1, 314.2) c1 = 7.304 (7.291, 7.317)
gof1 = struct with fields:
sse: 0.1770 rsquare: 0.9996 dfe: 898 adjrsquare: 0.9996 rmse: 0.0140
figure(4)
plot(Vcfit, tc, Vc), ylim([-1.5 1.5]), grid on
xlabel('Time (seconds)'), ylabel('V_{c}')
  3 Comments
Sam Chak
Sam Chak on 2 Nov 2023
Do you mean a single mathematical function that describes all 3 sine waves?

Sign in to comment.

More Answers (2)

Sam Chak
Sam Chak on 2 Nov 2023
I'm not an expert, but I have friends who are, and one of them said this is possible with this math function:
where
, and .
If you find the proposed math function useful, don't forget voting 👍 for the answer.
load('V_I.mat')
Vtime = Vabc3.Time;
Vdat1 = Vabc3.Data(:,1);
ta = Vtime(1:1001);
Va = Vdat1(1:1001);
tb = Vtime(1002:2100);
Vb = Vdat1(1002:2100);
tc = Vtime(2101:3001);
Vc = Vdat1(2101:3001);
subplot(211)
plot(Vtime, Vdat1, 'color', '#528AFA'), ylim([-1.5 1.5]), grid on
title('Plot from Data')
xlabel('Time (seconds)'), ylabel('V_{abc}')
[Vafit, gof1] = fit(ta, Va, 'sin1');
vaAmp = Vafit.a1; vaFreq = Vafit.b1; vaPhi = Vafit.c1;
[Vbfit, gof1] = fit(tb, Vb, 'sin1');
vbAmp = Vbfit.a1; vbFreq = Vbfit.b1; vbPhi = Vbfit.c1;
[Vcfit, gof1] = fit(tc, Vc, 'sin1');
vcAmp = Vcfit.a1; vcFreq = Vcfit.b1; vcPhi = Vcfit.c1;
t = linspace(0, 0.3, 3001);
t1 = 0.1;
t2 = 0.21;
Vaf = vaAmp*sin(vaFreq*t + vaPhi);
Vbf = vbAmp*sin(vbFreq*t + vbPhi);
Vcf = vcAmp*sin(vcFreq*t + vcPhi);
fun1 = sign(t - t1).*(Vaf/4 - Vbf/4); % sigma1
fun2 = Vaf/4 + Vbf/4 + Vcf/2 - fun1;
fun3 = Vaf/4 + Vbf/4 - Vcf/2 - fun1;
f = fun2 - sign(t - t2).*fun3; % <-- this math function
subplot(212)
plot(t, f, 'color', '#FA477A'), ylim([-1.5 1.5]), grid on
title('Plot using Math function')
xlabel('Time (seconds)'), ylabel('V')
  4 Comments
Atefeh
Atefeh on 13 Dec 2023
Dear Sam
Would you mind telling me the reference of the attached formulas?
Sam Chak
Sam Chak on 14 Dec 2023
The provided signum-based formula does describe the behavior of the piecewise function. However, I find it unnecessary to cite the formula, as it involves a clever but unintuitive manipulation of the sign functions to separate the sub-functions at different intervals in the domain.
Moreover, I have provided a more intuitively understandable formula in my Answer below. If you still wish to cite, consider referencing the special function used in the formulation, acknowledging the mathematical properties of that function.

Sign in to comment.


Sam Chak
Sam Chak on 14 Dec 2023
Given the piecewise function defined for
I propose a more intuitively understandable Heaviside-based formula, which I affectionately refer to as the Piecewise Function Put Together (PFPT) formula:
with .
In this formulation, the 1st term represents throughout the entire range of x. The 2nd term is activated by the Heaviside step function when x equals , plotting and canceling out from ​ onwards. Upon reaching , the Heaviside function triggers the 3rd term, plotting and canceling out from onwards.
If you appreciate the explanation of the PFPT formula, don't forget to vote 👍 for the answer as a token of appreciation.
load('V_I.mat')
Vtime = Vabc3.Time;
Vdat1 = Vabc3.Data(:,1);
ta = Vtime(1:1001);
Va = Vdat1(1:1001);
tb = Vtime(1002:2100);
Vb = Vdat1(1002:2100);
tc = Vtime(2101:3001);
Vc = Vdat1(2101:3001);
subplot(211)
plot(Vtime, Vdat1, 'color', '#528AFA'), ylim([-1.5 1.5]), grid on
title('Plot from Data')
xlabel('Time (seconds)'), ylabel('V_{abc}')
[Vafit, gof1] = fit(ta, Va, 'sin1');
vaAmp = Vafit.a1; vaFreq = Vafit.b1; vaPhi = Vafit.c1;
[Vbfit, gof1] = fit(tb, Vb, 'sin1');
vbAmp = Vbfit.a1; vbFreq = Vbfit.b1; vbPhi = Vbfit.c1;
[Vcfit, gof1] = fit(tc, Vc, 'sin1');
vcAmp = Vcfit.a1; vcFreq = Vcfit.b1; vcPhi = Vcfit.c1;
t = linspace(0, 0.3, 3001);
t1 = 0.1;
t2 = 0.21;
Vaf = vaAmp*sin(vaFreq*t + vaPhi); % f1, function at the 1st interval, t < t1
Vbf = vbAmp*sin(vbFreq*t + vbPhi); % f2, function at the 2nd interval, t1 < t < t2
Vcf = vcAmp*sin(vcFreq*t + vcPhi); % f3, function at the 3rd interval, t2 < t
%% Old Piecewise Function formula
% fun1 = sign(t - t1).*(Vaf/4 - Vbf/4); % sigma1
% fun2 = Vaf/4 + Vbf/4 + Vcf/2 - fun1;
% fun3 = Vaf/4 + Vbf/4 - Vcf/2 - fun1;
% fold = fun2 - sign(t - t2).*fun3; % <-- this math function
%% New Piecewise Function Put Together (PFPT) formula
fnew = Vaf + (Vbf - Vaf).*heaviside(t - t1) + (Vcf - Vbf).*heaviside(t - t2);
subplot(212)
plot(t, fnew, 'color', '#FA477A'), ylim([-1.5 1.5]), grid on
title('Plot using PFPT formula')
xlabel('Time (seconds)'), ylabel('V')
  1 Comment
Atefeh
Atefeh on 14 Dec 2023
Dear Sam
I really appreciate your time, it was more understandable.
Thanks a bunch

Sign in to comment.

Categories

Find more on Financial Toolbox in Help Center and File Exchange

Community Treasure Hunt

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

Start Hunting!