托伯利兹方程组解法
procedure TOEPLZ(R:array of double;var X:array of real;
Y:array of real; N:integer);
label 99;
var
G, H:array[0..100] of double;
M,M1,J,K,M2:integer;
SXN,SD,SGN,SGD,SHN,PT1,PT2,QT1,QT2,PP,QQ:double;
begin
If R[N] = 0 Then GoTo 99;
X[1]:=Y[1] / R[N];
If N = 1 Then Exit;
G[1]:=R[N - 1] / R[N];
H[1]:=R[N + 1] / R[N];
For M:=1 To N do
begin
M1:=M + 1;
SXN:=-Y[M1];
SD:=-R[N];
For J:=1 To M do
begin
SXN:=SXN + R[N + M1 - J] * X[J];
SD:=SD + R[N + M1 - J] * G[M - J + 1];
end;
If SD = 0 Then GoTo 99;
X[M1]:=SXN / SD;
For J:=1 To M do
X[J]:=X[J] - X[M1] * G[M - J + 1];
If M1 = N Then Exit;
SGN:=-R[N - M1];
SHN:=-R[N + M1];
SGD:=-R[N];
For J:=1 To M do
begin
SGN:=SGN + R[N + J - M1] * G[J];
SHN:=SHN + R[N + M1 - J] * H[J];
SGD:=SGD + R[N + J - M1] * H[M - J + 1];
end;
If (SD = 0) Or (SGD = 0) Then GoTo 99;
G[M1]:=SGN / SGD;
H[M1]:=SHN / SD;
K:=M;
M2:=(M + 1) div 2;
PP:=G[M1];
QQ:=H[M1];
For J:=1 To M2 do
begin
PT1:=G[J];
PT2:=G[K];
QT1:=H[J];
QT2:=H[K];
G[J]:=PT1 - PP * QT2;
G[K]:=PT2 - PP * QT1;
H[J]:=QT1 - QQ * PT2;
H[K]:=QT2 - QQ * PT1;
K:=K - 1;
end;
end;
ShowMessage('Never get here');
99: ShowMessage('Levinson method fails: singular principal minor');
end;