Bode Plot to Transfer Function

Question

Muhammad 2023-7-27

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/2001237-bode-plot-to-transfer-function

评论： Star Strider 2023-7-28

Hello, I take experimental data of my plant (Which is Atomic Force Microscopy PZT Actuator and Frame).

I have .DAT File, which contains following information Freuency (Hz), Magnitude dB, Phase Degree.

I want to estimate the transfer function of my AFM system from that bode plot, How can I do this in the MATLAB. I am not sure which order my plant is and either it is linear or non linear, I am new to the control system.

I have attached the .csv file of my data and Bode Plot image.

I

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

请先登录，再进行评论。

请先登录，再回答此问题。

Answer 1

Star Strider 2023-7-27

2
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/2001237-bode-plot-to-transfer-function#answer_1279597

在 MATLAB Online 中打开

If you have the System Identification Toolbox, this is (relatively) straightforward, however your data requires a bit of pre-processing. Start with the idfrd function and then use tfest. (The Signal Processing Toolbox has similar functions, however I usually use the System Identification Toolbox functions for these problems).

Try this —

figure

imshow(imread('Bode Plot_AFm.JPG'))

T1 = readtable('HH8.csv')

T1 = 1001×4 table

F GdB P G ___ ____ ____ ____ 100 33.9 -110 49.4 101 34.1 -110 50.7 101 33.8 -107 48.9 102 34.2 -108 51.1 102 33.6 -105 47.8 103 34 -107 49.8 103 33.1 -111 45.2 104 33.8 -108 48.8 104 33.7 -108 48.4 105 33.5 -106 47.5 105 33.7 -106 48.5 106 33.6 -107 48 106 34 -105 50.1 107 34 -103 50 107 34.1 -103 51 108 34.2 -102 51.6

Ts = 0; % Fill With Actual Sampling Frequency

FHz = T1.F;

for k = 1:size(T1,1)-1

if FHz(k+1) == FHz(k)

FHz(k+1) = FHz(k+1)+0.5; % 'Brute Force' Interpolation

end

Mag = T1.G;

PhDeg = T1.P;

Response = Mag.*exp(1j*deg2rad(PhDeg)); % Complex Vector

sysfr = idfrd(Response, FHz, Ts, 'FrequencyUnit','Hz')

sysfr = IDFRD model. Contains Frequency Response Data for 1 output(s) and 1 input(s). Response data is available at 1001 frequency points, ranging from 100 Hz to 1.5e+04 Hz. Status: Created by direct construction or transformation. Not estimated.

tfsys = tfest(sysfr,18)

tfsys = -6.838e04 s^17 + 1.058e10 s^16 - 5.832e14 s^15 + 1.29e19 s^14 - 2.521e23 s^13 + 5.075e27 s^12 - 2.423e31 s^11 + 9.36e34 s^10 - 9.093e37 s^9 + 3.076e41 s^8 - 1.342e44 s^7 + 4.287e47 s^6 - 9.57e49 s^5 + 2.978e53 s^4 - 3.277e55 s^3 + 1.019e59 s^2 - 4.296e60 s + 1.372e64 -------------------------------------------------------------------------------------------------------------------------------------------------------------------------- s^18 + 2.677e05 s^17 + 4.519e09 s^16 + 3.152e14 s^15 + 1.125e18 s^14 + 8.793e22 s^13 - 3.708e26 s^12 + 2.531e30 s^11 - 2.069e33 s^10 + 8.788e36 s^9 - 4.455e39 s^8 + 1.255e43 s^7 - 4.823e45 s^6 + 8.854e48 s^5 - 2.802e51 s^4 + 3.061e54 s^3 - 8.38e56 s^2 + 4.148e59 s - 1.02e62 Continuous-time identified transfer function. Parameterization: Number of poles: 18 Number of zeros: 17 Number of free coefficients: 36 Use "tfdata", "getpvec", "getcov" for parameters and their uncertainties. Status: Estimated using TFEST on frequency response data "sysfr". Fit to estimation data: 85.98% FPE: 4.104, MSE: 3.819

figure

compare(sysfr, tfsys)

Experiment to get desired results. Having the actual sampling frequency would likely improve this.

.

6 个评论
显示 4更早的评论隐藏 4更早的评论

Star Strider 2023-7-28

That depends on the data you have. The easiest way would be to check the device sampling rate, assuming that it is constant. If you have time-domain data, and have the associated time vector, checking the sampling rate is easier, since this involves taking the inverse of the mean of the sampling time differences.

That would go something like this:

Ts = mean(diff(t))

Fs = round(1/Ts)

[sr,tr] = resample(s, t, Fs);

assuming here that ‘t’ is the time vector and ‘s’ is the signal vector.

The resample call is necessary if the standard deviation of the time differences is greater than about

of ‘Ts’ since that would indicate significant non-uniform sampling intervals, and all signal processing (except the nufft function) require that the signal be uniformly sampled. So that is relatively important. You would then use the resampled time-domain vectors with iddata (instead of idfrd). The rest of the code remains essentially unchanged.

Having the sampling interval ‘Ts’ and knowing that it is uniform across the signal, would likely make the system identification easier and more reliable.

One other way to estimate the sampling interval (again assuming that it is constant) and assuming the highest frequency in the observed data in the Bode plot is the Nyquist frequency (one-half the sampling frequency) would be to simply take the inverse of the Nyquist frequency and divide that result by 2. (It would be necessary to be certain that all these assumptions are true first. Otherwise, the resulting sampling interval calcuation would be completely unreliable.)

.

Star Strider 2023-7-28

在 MATLAB Online 中打开

HH8.csv

A leap-of-faith is involved in assuming that the highest frequency is the Nyquist frequency, however using that, the result is reasonable if not perfect —

T1 = readtable('HH8.csv')

T1 = 1001×4 table

F GdB P G ___ ____ ____ ____ 100 33.9 -110 49.4 101 34.1 -110 50.7 101 33.8 -107 48.9 102 34.2 -108 51.1 102 33.6 -105 47.8 103 34 -107 49.8 103 33.1 -111 45.2 104 33.8 -108 48.8 104 33.7 -108 48.4 105 33.5 -106 47.5 105 33.7 -106 48.5 106 33.6 -107 48 106 34 -105 50.1 107 34 -103 50 107 34.1 -103 51 108 34.2 -102 51.6

FHz = T1.F;

Ts = 1/(2*max(FHz)) % Fill With Actual Sampling Frequency

Ts = 3.3333e-05

for k = 1:size(T1,1)-1

if FHz(k+1) == FHz(k)

FHz(k+1) = FHz(k+1)+0.5; % 'Brute Force' Interpolation

end

Mag = T1.G;

PhDeg = T1.P;

Response = Mag.*exp(1j*deg2rad(PhDeg)); % Complex Vector

sysfr = idfrd(Response, FHz, Ts, 'FrequencyUnit','Hz')

sysfr = IDFRD model. Contains Frequency Response Data for 1 output(s) and 1 input(s). Response data is available at 1001 frequency points, ranging from 100 Hz to 1.5e+04 Hz. Sample time: 3.3333e-05 seconds Status: Created by direct construction or transformation. Not estimated.

tfsys = tfest(sysfr,5) % Experiment With The Pole and Zero Arguments

Warning: The frequency response of the data is not symmetric. The estimated model may have complex coefficients.

tfsys = (-1.074e04 + 0.71i) s^4 + (-2.52e09 - 9.5e02i) s^3 + (3.935e12 + 2.7e08i) s^2 + (-57049681860136032 - 8.5e11i) s + (4.185e21 + 4e07i) ---------------------------------------------------------------------------------------------------------------------------------------------- s^5 + (3.181e04 - 0.074i) s^4 + (9.282e08 + 2e02i) s^3 + (1.361e13 - 3.2e07i) s^2 + (1.153e17 - 1.3e09i) s + (-27042577579097112576 + 5.8e09i) Continuous-time identified transfer function. Parameterization: Number of poles: 5 Number of zeros: 4 Number of free coefficients: 10 Use "tfdata", "getpvec", "getcov" for parameters and their uncertainties. Status: Estimated using TFEST on frequency response data (complex) "sysfr". Fit to estimation data: 77.47% FPE: 10.05, MSE: 9.854

figure

compare(sysfr, tfsys)

This will require some fine-tuning, however it will get you started. The AFM actuator is likely a cantilever, however I do not know what sort of dynamics the frame could have. The tfest function allows separately sepcifying the numbers of zeros as well as the numbers of poles, so experiment with that as well. (It might be easier to initially model this in state-space (use ssest) and then convert it to a transfer function. I leave that to you to explore.)

.

Muhammad 2023-7-28

移动：Star Strider 2023-7-28

Thank You so much.

Star Strider 2023-7-28

As always, my pleasure!

请先登录，再进行评论。

Bode Plot to Transfer Function

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

采纳的回答

6 个评论
显示 4更早的评论隐藏 4更早的评论

更多回答（0 个）

另请参阅

类别

标签

产品

版本

Community Treasure Hunt

Bode Plot to Transfer Function

0 个评论 显示 -2更早的评论隐藏 -2更早的评论

采纳的回答

6 个评论 显示 4更早的评论隐藏 4更早的评论

更多回答（0 个）

另请参阅

类别

标签

产品

版本

Community Treasure Hunt

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

6 个评论
显示 4更早的评论隐藏 4更早的评论