How to make a transfer function minimum phase?

Dear MATLAB Community,

I have a plant Transfer Function which is non minimum phase. I want to make it stable minimum phase system so that I can inverse it without instaability.

% Define the transfer function

H = tf([-4.8 16000 0 0],[4.8 16080 286800 51160000]);

isminphase([-4.8 16000 0 0], [4.8 16080 286800 51160000])

Thanks!

[Paul]

Hi Govind,

H is defined as a continous time transfer function (Control System Toolbox), but isminphase as used above is a Signal Processing Toolbox function that is only applicable to discrete time filters.

Assuming that H really is a continous time transfer function, it isn't minimum phase. But there is no recipe "to make it stable, minimum phase." Such a goal is too open ended. There would have to be some other criteria as well.

[Sam Chak]

Hi @Govind Narayan Sahu

What do you want to do after getting the inverse of this function?

H = tf([-4.8 16000 0 0], [4.8 16080 286800 51160000])

H =
 
            -4.8 s^3 + 16000 s^2
  -----------------------------------------
  4.8 s^3 + 16080 s^2 + 286800 s + 5.116e07
 
Continuous-time transfer function.

step(H, 0.01)

[Govind Narayan Sahu]

Hi @Sam Chak,

I want to design a transfer function (Hinv) that will behave inversely to H so that the final magnitude of H*Hinv is (approximately) equal to one. This is because when I pass an input signal, it amplifies the output at about 10 Hz.

H will be invertible if it has a minimum phase.

I read this article: https://doi.org/10.1080/10910344.2012.698970

I am able to create a frequency response that has a minimum phase, but not able to get a correct transfer function with the use of "invfreqz" and "freqz".

% Define the transfer function

num = [-4.8 16000 0 0];

den = [4.8 1.608e04 2.868e05 5.116e07];

H = tf(num,den); % Inverse Transfer function Den/Num

w = 2*pi*(0.1:0.01:100);

% Calculate complex cepstrum

c_hat_n = ifft(log(abs(squeeze(freqresp(H,w)))));

N = length(c_hat_n);

Nst = round(N/2)+1;

c_hat_n(Nst:N) = 0;

c_hat_n = c_hat_n*2;

h_hat_min = c_hat_n;

Hmin=exp(fft(h_hat_min));

w11 = angle(Hmin);

[b,a] = invfreqz((Hmin),w11, 3,3, [], 1000);

[h,w1] = freqz(b,a,N);

subplot(211)

hold on;plot(w/2/pi, abs(squeeze(freqresp(H,w))));

plot(w/2/pi, abs(Hmin),'r')

plot(w/2/pi, abs(h),'g')

xlabel('Frequency [Hz]')

ylabel('Magnitude, [N/N]')

subplot(212)

hold on;plot(w/2/pi, phase(squeeze(freqresp(H,w)))+2*pi);

plot(w/2/pi, phase(Hmin),'r')

plot(w/2/pi, 2*pi-w1,'g')

xlabel('Frequency [Hz]')

ylabel('Phase, [degree]')