From time to time, a control engineer faces a problem in opposite direction. For example, she obtains a discrete-time plant from identification process and wants to convert it to continuous-time domain for analysis and design. Surprisingly enough, I could not find a Scilab command for such purpose. Maybe … I don’t know, it might appear in future version. But at this time I have to get around this problem by writing my own dls2cls function*. Bilinear transform method is used. * Actually I prefer a shorter command name like d2c, but well, dls2cls makes it recognized as an exact inverse of cls2dls. Before coding anything, I google to find if there exists a solution already. Could not find one other than some suggestion to use horner command in Scilab to convert polynomial of s to z using the Bilinear transform relationship. Conceptually, that should work. However, I tested it with a stable 5-th order transfer function and found that the result came out having higher order and, worse, an unstable continuous-time system. So this approach is not reliable due to some Scilab internal computational error beyond my comprehension. I have no choice but to write my own algorithm from scratch. So here it goes. Start from the Bilinear transform relationship
(1)where is sampling period. From (1), we derive for in terms of as
(2)To make the notation concise, define
(4)From now on we simply write them as and . So the relationship (2) becomes
(5)To demonstrate the process, let us consider a third order discrete-time transfer function of the form
(6)Then its continuous-time counterpart can be found using (5) as
(7)Multiplying both numerator and denominator by yields
(8)This derivation can be extended to an order discrete-time transfer function
(9)One can show that the corresponding continuous-time transfer function is
(10)The equation (10) can be easily written as Scilab algorithm. Indeed, this is a good exercise on writing a Scilab function. We show it step by step below. Create a new file in Scinotes. Start with Scilab function placeholder.
Save the file as dls2cls.sci in your chosen directory. The algorithm will sit between the function and endfunction keyword. The arguments passed to dls2cls are a discrete-time transfer function in domain and sampling time . The function returns a continuous-time transfer function Sysc. At the beginning of function, define and variable
function Sysc=dls2cls(Sysd, T) endfunction
While dls2cls should work with any valid discrete-time system, we restrict it only to a proper (which includes strictly-proper) rational function; i.e., the denominator order must be greater or equal to the numerator order. This represents a system that is realizable. Use the following code to perform this check.
The error command displays the quoted message and exits the function in case an improper transfer function is passed to dls2cls. To compute the algorithm (10), define and
numdeg = degree(Sysd.num); dendeg = degree(Sysd.den); if numdeg>dendeg error('The discrete-time system is improper'); end
and extract the coefficients from numerator and denominator
f = -(s*T+2); g = (s*T-2);
Now, if the transfer function is strictly proper, the length of bn is shorter than an. Then we have to pad bn with 0 so that they have equal length.
bn = coeff(Sysd.num); an = coeff(Sysd.den);
The essence of this algorithm is the for loop that implements (10)
npad = dendeg - numdeg; // number of zero to be padded, if any if npad > 0 for i=0:npad bn = [bn 0]; // pad numerator coeff vector with zero end end
The computation is straightforward. After we derive the numerator and denominator of (10), a continuous-time transfer function is formed using syslin command
nums = 0; dens = 0; for i=0:dendeg, nums = nums + bn(i+1)*(f^i)*(g^(dendeg-i)); dens = dens + an(i+1)*(f^i)*(g^(dendeg-i)); end
and returned to the caller.
Sysc = syslin('c',nums,dens);
Download the complete file –> dls2cls.sci
To test the function, load it into Scilab workspace
then create an arbitrary discrete-time model, say
Then pass to dls2cls with some specified sampling time , say, 0.01
-->numd = -2*z^4-z^3+z^2+2*z+3; -->dend = 3*z^4+2*z^3+z^2 -z -2; -->Sysd = syslin('d',numd/dend) Sysd = 2 3 4 3 + 2z + z - z - 2z ------------------- 2 3 4 - 2 - z + z + 2z + 3z
The function returns the continuous-time conversion. To verify the result, convert it back to discrete-time using Scilab cls2dls command. Note that this command requires a system with state-space representation.
-->Sysc=dls2cls(Sysd,0.01) Sysc = 2 3 4 48 - 2.08s + 0.0016s - 0.000028s + 1.000D-08s --------------------------------------------- 2 3 4 48 + 2.08s + 0.0016s + 0.000028s + 1.000D-08s
Then converts back to transfer function
-->Syscss=tf2ss(Sysc); -->Sysdss = cls2dls(Syscss,0.01);
We get the same discrete-time transfer function back (only normalized by dividing the numerator and denominator by 3). I am happy with the result, though I have not yet tested this function extensively. If you find a bug kindly let me know.
-->Sysd1 = ss2tf(Sysdss) Sysd1 = 2 3 4 1 + 0.6666667z + 0.3333333z - 0.3333333z - 0.6666667z ----------------------------------------------------- 2 3 4 - 0.6666667 - 0.3333333z + 0.3333333z + 0.6666667z + z
Download the complete file –> dls2cls.sci