matlab system identification toolbox gives me incorrect DC gain of transfer function

Question

everyone. I'm trying to build a noise model using system identification toolbox.

My system can be written as: y = C(z)/D(z) e . Basically I measure system response y purely due to unmeasured white noise e . I think system id will do the job since it's just an ARMA model.

Before applying it to real data, I wrote the following test script, which simulates a known system with whitenoise input and try to use the response y to get an estimated model and compare it with the known one by plotting their bode. My problem is that the estimated bode has the same shape as the true one, but the DC gain is very different. Can some one read my script and tell me what is wrong? THANK YOU!

close all; close all;
wn = 10;
sys = tf([wn^2], [1, 2*0.9*wn wn^2]);

% discretize;
Ts    = 0.025;
sysdt = c2d(sys, Ts);
Fs    = 1/Ts;
Tn    = 4;
N     = Tn/Ts;
Tsim = (0:N-1)' * Ts;

whitenoise = randn(N, 1);
Y          = lsim(sysdt, whitenoise, Tsim);

%the "input" u is 0, we have only noise e
td_data    = iddata(Y, zeros(size(Y)), Ts);

%% estimate use different commands
%1.armax model: A(q) y(t) = B(q) u(t-nk) +  C(q) e(t)
% syntax: M = armax(data,  [na, nb, nc, nk]
idout_armax = armax(td_data, [2, 0, 1, 1]);
idsys_armax = tf(idout_armax.c, idout_armax.a, Ts);
figure, bode(sysdt, idsys_armax)
legend('true model', 'armax')

%2. Box_Jenkins:  y(t) = [B(q)/F(q)] u(t-nk) +  [C(q)/D(q)]e(t)
%                      b+1 c  d  f  k
idout_bj = bj(td_data, [1, 1, 2, 1, 0]);
idsys_bj = tf(idout_bj.c, idout_bj.d, Ts);
figure, bode(sysdt, idsys_bj)
legend('true model', 'box jenkins')

%3. If I use the whitenoise data as input *u* , I can get correct DC gain with oe (most of the time).
td_data_wn = iddata(Y, whitenoise, Ts);
% OE model: y(t) = [B(q)/F(q)] u(t-nk) +  e(t)
%                         nb nf nk
idout_oe = oe(td_data_wn, [1, 2, 0]);
idsys_oe = tf(idout_oe.b, idout_oe.f, Ts);
figure, bode(sysdt, idsys_oe), legend('sysdt', 'idsys oe')

Answer 1

I found out a possible reason for the DC gain by myself:

The system identification toolbox gives an estimation of the noise variance. So although I generate my data using randn with variance 1, the toolbox assumes (estimates) another noise variance, which I should use to scale my estimated transfer function to get the right DC gain.

So if in the above code, the armax model estimation part, I use

idout_armax = armax(td_data, [2, 0, 1, 1]);
scale = sqrt(idout_armax.NoiseVariance);
idsys_armax = tf(idout_armax.c, idout_armax.a, Ts)*scale;
figure, bode(sysdt, idsys_armax)

the DC gain should (almost) match sometimes.

I hope this is the correct reason.