Matlab 测试程序如下:
Ts = 0.05;
sysc = tf([0 1], [1 1]);
sysd = c2d(sysc, Ts);
t = 0.1:Ts:20;
x = square(t);
yd = lsim(sysd, x, t, 0);
ys = yd + 0.1*rand(size(yd));
[numc, denc] = LeastSquare(x, ys, 60);
sysd_hat = tf(numc, denc, Ts);
yhat = lsim(sysd_hat, x, t, 0);
plot(t, ys, t, yhat);
最小二乘法程序, x, y 为输入输出, N为阶数:
function [numc, denc] = LeastSquare(x, y, N)
count = length(y);
M = count - 1;
ai = zeros(N*2, M);
for i=1:N
ai(i, i:M) = x(1:(count-i));
ai(i+N, i:M) = -y(1:(count-i));
end
bi = y(2:(count));
xd = (inv(ai*ai')*ai*bi)';
numc = [0 xd(1:N)];
denc = [1 xd(N+1: N+N)];
|