From bode to transfer function

I'm trying to get a transfer function out of the bode plot data. I've tried multiple functions(using this and this) yet when I plot both an example function and the solution the functions give me, it doesn't match.

So first I made the example:

examp1=tf([1],[1 -4 3])
[AMP PHA W]=bode(examp1)

which is a function with 2 poles. Then with tfest I made(copy/paste):

gain = squeeze(AMP);
phase = squeeze(PHA);
whz = squeeze(W);

response = gain.*exp(1i*phase*pi/180);
Ts = 0.1; % your sampling time
w=whz*2*pi; %convert Hz to rad/sec 
gfr = idfrd(response,w,Ts);

sys=tfest(gfr,2);
bode(sys)
hold on
bode(examp1)

within the tfest function the 2nd value should be 2, correct? Since we have 2 poles in the original function. Yet in the end, when we compare the 2 bodes, they aren't much alike.

How should I do this?

Answer 1

[Arkadiy Turevskiy]

You need to set tfestoptions Focus property to "prediction" to allow for unstable systems (your system has two unstable poles).

Doc page explaining that is here, look under "Focus" property.

Here is the code that gives a perfect fit (also note that bode returns frequency vector in rad/time unit, so there is no need to convert to hz as you were doing).

examp1=tf([1],[1 -4 3]);
[AMP PHA W]=bode(examp1);

gain = squeeze(AMP);
phase = squeeze(PHA);
w = squeeze(W);

response = gain.*exp(i*phase*pi/180);
Ts = 0.1; % your sampling time 
gfr = idfrd(response,w,Ts);

Options = tfestOptions;                             
Options.Focus = 'prediction';   

tf1 = tfest(gfr, 2, 1, Options, 'Ts', 0.1);
P = bodeoptions;
P.PhaseWrapping = 'on';

bodeplot(examp1,tf1,P);

[Raymond]

Thank you. 'gfr' is indeed a perfect fit. If I create a bodeplot of 'tf1', (which is the predicted transfer function, correct?) the values slightly rise near the end, untill theres a black line in the plot. Is this because of the unstable poles?

In the future we plan to do experiments with an yet to determine system. It is uncertain how many poles and zero's the system has. How should we use the options of tfest then?

Thanks again.

[Raymond]

By the way, why does the system have 2 unstable poles? Aren't the poles just 3 and 1, in which case they are stable?

[Arkadiy Turevskiy]

I meant to compare tf1 with examp1 in my response, but typed in the wrong variable name, thanks for noticing this. I corrected the answer, and the fit is still almost perfect. The reason why it is not completely dead on is because you are trying to estimate a discrete-time transfer function. If you estimate a continuous one, the fit is even better:

examp1=tf([1],[1 -4 3]);
[AMP PHA W]=bode(examp1);
gain = squeeze(AMP);
phase = squeeze(PHA);
w = squeeze(W);
response = gain.*exp(i*phase*pi/180);
sys=frd(response,w);
gfr=idfrd(sys);
Options = tfestOptions;
Options.Focus = 'prediction';
tf2=tfest(gfr,2,1,Options);
P = bodeoptions;
P.PhaseWrapping = 'on';
bodeplot(examp1,tf2,P)

[Arkadiy Turevskiy]

As for the other question, about now knowing the number of poles and zeros beforehand, there is no "sceintific" way to come up with best options - try a bunch of different options and number of poles and zeros and pick the one combination that gives the best fit. You can use compare to automate this.

You may also look at estimating a state space instead of a transfer function, using ssest. When you call this function, you can specify system order as a vector, say [1 10], and the function will then return a plot helping you choose the best order as shown here .

If needed, you can then convert the identified state-s[ace model into a transfer function using tf .

[Arkadiy Turevskiy]

Finally, the system has two unstable poles because the poles (at 1 and 3) are both positive numbers. The pole is stable if real part is less than zero, and unstable if it is >0, which is your case. I am not going to explain why- read wikipedia or take a look at any introductory controls textbook.

To convince yourself in MATLAB quickly:

s=tf('s');
sys1=1/(s-1); %positive pole at 1
step(sys1); %unstable response

Now for the negative pole:

sys2=1/(s+1); % pole at -1
close all;
step(sys2); %stable response