Bode Plot to Transfer Function

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

Answer 1

[Star Strider]

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
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.

.

[Muhammad]

Thanks @Star Strider for providing me code with all steps, I am happy to see this code is working for me. One thing I just want to know, how to choose my sampling frequency from my data or I have to check my device sampling rate? I know it's not the right platform to ask this question.

[Star Strider]

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.)

.

[Muhammad]

@Star Strider My another question is How to choose sampling rate based on my bode plot? and why you wrote 18 in the "tfsys = tfest(sysfr,18)"

when I change sampling rate to some other value like

Ts = 0.01; or Ts = 0.1; I get following error.

Error using tfest (line 113)

Index exceeds the number of array elements (0).

Error in bodetf (line 16)

tfsys = tfest(sysfr,18)

and when I change it to 0.0001 the mode appears to be just -7 correct and when I further decrease it to 0.0001 mode firs only 30%....

[Star Strider]

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
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]

Thank You so much.

[Star Strider]

As always, my pleasure!