Effacer les filtres
Effacer les filtres

discrete to continuous conversion of transfer function

15 vues (au cours des 30 derniers jours)
Pierce Hart
Pierce Hart le 31 Août 2017
Hello, I have the following discrete transfer function given in the first image. the values for each term are as given:
  • a1= -2.9964e+00
  • b1= -9.2542e-03
  • a2= 2.9933e+00
  • b2= 1.9492e-02
  • a3= -9.9694e-01
  • b3= -1.0238e-02
Using the old thd2thc command the continuous time transfer function should be as the same as shown in figure 2. The sampling time is 2e-8
My code to tackle this problem is as follows. The last term on the numerator and last 2 on the denominator are out by large values.
clear all;
clc;
numdsp=[0,-9.2542e-3,1.9492e-2,-1.0238e-2];
dendsp=[1,-2.9964,2.9933,-9.9694e-1];
Hdsp=filt(numdsp,dendsp,2e-8)
Hc=d2c(Hdsp)
Is the difference due to the old thd2thc method or is there an issue with where I'm going.

Réponses (1)

Emma Farnan
Emma Farnan le 8 Nov 2022
Unfortunately, I expect that this is too late for the response to be useful for you but I did run into the same issue when investigating the 1996 paper by Rotea et al. that these values come from.
It is my understanding that the thd2thc command is significantly outdated, and has been replaced with d2c for discrete to continuous transfer function conversions. That being said, I don't think that there is any issue with the replacement of the old to new version of this MATLAB function and is simply a precision error caused by the number of decimals included in a1-b3 that you posted. Having 5 significant figures sounds sufficient, but because some of these poles and zero locations are so close to the instability break points on the pole-zero maps, the rounding that you used causes the transfer function to just barely fall into the unstable region. You can prove this by taking Rotea's results after the discrete-to-continous conversion and working backwards to obtain discrete values for comparison...
% NOTATION: ---------------------------------------------------------
% H_ = Transfer Functions originating with the b/a coefficients
% h_ = Computed TF starting with the published Continous results
% _d = Discrete-Time TF
% _c = Continous-Time TF
% _zpk = Zero-pole-Gain TF Format (similar to form in Eq 16 and 17 of
% Rotea et al. paper)
% -------------------------------------------------------------------
% Set the sampling frequency and period for the discrete TF
fs = 50e6; % 50 MHz Sampling Frequency
Ts = 1/fs; % associated sampling period
% Input the coefficients found for equation (11) in Rotea et al.
b1 = -9.2542e-3; b2 = 1.9492e-2; b3 = -1.0238e-2; % Numerators
a1 = -2.9964e0; a2 = 2.9933e0; a3 = -9.9694e-1; % Denominators
% Format num/den form with k = 0 terms
num = [0 b1 b2 b3];
den = [1 a1 a2 a3];
% Generate the discrete TF in digital signal processing (DSP) form
% This would be the same as "tf(num,den,'Variable','z^-1')"
Hd = filt(num,den,Ts);
% Convert this discrete-time transfer function to Continuous-time TF
Hc = d2c(Hd);
% Also convert to zero-pole-gain format
[z,p,k] = tf2zp(Hc.num{1},Hc.den{1});
Hzpk = zpk(z,p,k);
% Plot the step response and bode plots
figure('Name','Step Response from Rotea b/a coefficients')
step(Hc,'r',Hd,'b--',Hzpk,'g-.')
figure('Name','Bode Plots from Rotea b/a coefficients')
bode(Hc,'r',Hd,'b--',Hzpk,'g-.')
% -------------------------------------------------------------------------
% The above results appear to be unstable so work back from the published
% continous-time transfer function presented in Eq 16 of Rotea. It is a
% continuous-time TF presented in zero-pole-gain format
K = -4.8810e5; % Gain
Z = [-1.7152e4; 5.0618e6]; % Zeros of Continuous TF
P = [-3.2607e4; -6.0270e4-1.17e6j; -6.0270e4+1.17e6j]; % Poles of TF
% Generate zero-pole-gain TF model and work backwards
hzpk = zpk(Z,P,K);
% Convert to continuous-time TF
[bnum,aden] = zp2tf(Z,P,K);
hc = tf(bnum,aden);
% Finally Convert to Discrete-time TF
hd = c2d(hc,Ts);
% Plot these with a step response and bode plots as well
figure('Name','Step Response after working backwards')
step(hc,'r',hd,'b--',hzpk,'g-.')
figure('Name','Bode Plots after working backwards')
bode(hc,'r',hd,'b--',hzpk,'g-.')
% -------------------------------------------------------------------------
% Plot the pole-zero maps of the continous and discrete TFs
figure('Name','Pole-Zero Map of Continuous-Time TFs')
pzmap(Hc,'r',hc,'b')
legend('TF using d2c','Published TF after d2c')
figure('Name','Pole-Zero Map of Discrete-Time TFs')
pzmap(Hd,'r',hd,'b')
xlim([0.98 1.03]) % Zoom in to show how close to stable original is
legend('Original Discrete-Time TF','TF working backwards using c2d')
% -------------------------------------------------------------------------
% Print out the transfer functions
brkline = '==============================================================';
abs = '***** TF from plugging in b and a coefficients into eq (11)';
rev = '***** TF working backwards from zero-pole-gain TF in eq (16)';
% Discrete TFs
fprintf('%s\nDISCRETE-TIME TRANSFER FUNCTIONS\n%s\n',brkline,brkline)
============================================================== DISCRETE-TIME TRANSFER FUNCTIONS ==============================================================
fprintf('%s%s%s%s',abs,evalc('Hd'),rev,evalc('hd'))
***** TF from plugging in b and a coefficients into eq (11) Hd = -0.009254 z^-1 + 0.01949 z^-2 - 0.01024 z^-3 -------------------------------------------- 1 - 2.996 z^-1 + 2.993 z^-2 - 0.9969 z^-3 Sample time: 2e-08 seconds Discrete-time transfer function. ***** TF working backwards from zero-pole-gain TF in eq (16) hd = -0.009254 z^2 + 0.01949 z - 0.01024 ----------------------------------- z^3 - 2.996 z^2 + 2.993 z - 0.9969 Sample time: 2e-08 seconds Discrete-time transfer function.
% Continuous-time TFs
fprintf('\n%s\nCONTINUOUS-TIME TRANSFER FUNCTIONS\n%s\n',brkline,brkline)
============================================================== CONTINUOUS-TIME TRANSFER FUNCTIONS ==============================================================
fprintf('%s%s%s%s',abs,evalc('Hc'),rev,evalc('hc'))
***** TF from plugging in b and a coefficients into eq (11) Hc = -4.881e05 s^2 + 2.463e12 s - 2.504e16 ------------------------------------------ s^3 + 1.532e05 s^2 + 1.402e12 s - 5.008e18 Continuous-time transfer function. ***** TF working backwards from zero-pole-gain TF in eq (16) hc = -488100 s^2 + 2.462e12 s + 4.238e16 ---------------------------------------- s^3 + 153147 s^2 + 1.376e12 s + 4.475e16 Continuous-time transfer function.
% Zero-Pole-Gain TFs
fprintf('\n%s\nZERO-POLE-GAIN TRANSFER FUNCTIONS\n%s\n',brkline,brkline)
============================================================== ZERO-POLE-GAIN TRANSFER FUNCTIONS ==============================================================
fprintf('%s%s%s%s\n',abs,evalc('Hzpk'),rev,evalc('hzpk'))
***** TF from plugging in b and a coefficients into eq (11) Hzpk = -4.8811e05 (s-5.037e06) (s-1.018e04) --------------------------------------- (s-1.4e06) (s^2 + 1.553e06s + 3.577e12) Continuous-time zero/pole/gain model. ***** TF working backwards from zero-pole-gain TF in eq (16) hzpk = -4.881e05 (s+1.715e04) (s-5.062e06) ----------------------------------------- (s+3.261e04) (s^2 + 1.205e05s + 1.373e12) Continuous-time zero/pole/gain model.
% =======================================================================
% What seems peculiar at this point is that the printed versions for the
% continuous-time TF and the zpk-model TF differ whether we worked forwards
% or backwards, but the discrete-time TF appears to have the same
% coefficients whether we moved forwards or backwards. But investigating a
% little further, it can be seen that this is just due to the rounding on
% the values and this small rounding is what makes all of the difference
% Print the values of in the default rounded format
% First the numerators (top: Forward method. Bottom: Backwards method)
disp([num; hd.num{1}])
0 -0.0093 0.0195 -0.0102 0 -0.0093 0.0195 -0.0102
% Then the Denominators
disp([den; hd.den{1}])
1.0000 -2.9964 2.9933 -0.9969 1.0000 -2.9964 2.9933 -0.9969
% The two methods appear to match, but lets print more decimals
format long
% Again the numerator first
disp([num; hd.num{1}])
0 -0.009254200000000 0.019492000000000 -0.010238000000000 0 -0.009254180469808 0.019492064920648 -0.010237545970321
% Then the Denominators
disp([den; hd.den{1}])
1.000000000000000 -2.996400000000000 2.993300000000000 -0.996940000000000 1.000000000000000 -2.996391849588389 2.993333953072731 -0.996941746015171
So to sum it up... there is nothing wrong with d2c in comparison to thd2thc but for this specific example, small artifacts from rounding precision makes the original Discrete-time TF unstable and therefore the discrete-to-continuous conversion gives results that don't match what is expected.

Catégories

En savoir plus sur Get Started with Control System Toolbox dans Help Center et File Exchange

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by