关于SVD-TLS算法的总结

1. 必须重新列出法方程组，使它只包含p个未知数；
2. 求解Ax=b的最小二乘法只认为b含有误差，当然实际上矩阵A也是存在误差。
   故一种比最小二乘法更合理的方法是同时考虑A和b的误差和扰动。
   令m*n矩阵A的误差为E，向量b的误差为e，则矩阵方程的最小二乘解为：
   
   因为考虑了两者的误差与扰动故可称为总体最小二乘法。

我们直接上SVD-TLS算法的具体步骤：

步骤1：计算增广矩阵B的SVD，并存储奇异值和矩阵V；
步骤2：确定增广矩阵B的有效秩p；
步骤3：利用式（3.4.56）和式（3.4.53）计算矩阵S（p）；
步骤4：求S（p）的逆矩阵S-§，并由式（3.4.59）计算未知参数的总体最小二乘估计。

而对应的MATLAB函数可以为：

B=[-r',R];                                  %对应矩阵B
[U,K,V]=svd(B);              %求增广矩阵B的SVD，并存储奇异值和矩阵V
P=rank(B);                  % 求得增广矩阵B的有效秩，定义为P
S=zeros(P+1);          %构造一个(p+1)*(p+1)维的矩阵S                   

for j=1:p                  
for i=1:p+1-P
  djj=K(j,j)*K(j,j);     %构造djj,并求其平方
  vij=V(i:i+p,j);         %构造vij
 S= S+djj*vij*vij';     %计算矩阵S的二重级数求和
 end
end

Sni=inv(S);                   %求S逆矩阵
a=zeros(1,P);                     %构造a矩阵
for i=1:P                                
 a(1,i)=Sni(i+1,1)/Sni(1,1);     %求出矩阵AR模型参数a=[a1,...,ap]'
end

我们也可以带入时间序列确定其自相关函数，构造一个矩阵来运用SVD-TLS算法来求其功率谱密度：
设
这里可以看到两个余弦函数确定的时间序列给定了幅值和频率，其中V（n）为噪声函数，它可以由均值和方差的函数来确定。

clear all;                       %清除workspace之前的变量
close all;                       %关闭之前的图像
clc;                           %清除命令行之前的文字
n=[1:128];                     %取定采样点n=1至128
f1=0.2;                        %取定f1频率的值
f2=0.213;                      %取定f2频率的值
A=sqrt(20);                     %取定第一个正弦函数的振幅
B=sqrt(2);                      %取定第二个正弦函数的振幅
x=A*cos(2*pi*f1*n)+B*cos(2*pi*f2*n);     %定义x函数  
noise=0+1*randn(1,length(n));        %添加均值为0、方差为1的高斯白噪声
xn=x+noise;                %在x1基础上添加加性高斯白噪声，定义xn函数 
m=xcorr(xn);                    %m为xn的自相关函数（序列)

%-----------------------------------------

p=4;                    %取定R的阶数,更改p=4，64，100，的值，观察%相对应的谱估计
q=125;                    %此处一定要满足q>=p

for i=1:p
    for j=1:p
      R(i,j)=m(q+i+j-1-p);  %构造一个pxp阶的自相关矩阵（Hankel矩阵）
%(课本P88 3.4.33a)
 end
end
Rlegnth=size(R)                      %输出验证R矩阵的行列数的值


for i=1:p                                   %i=1~p
    r(i)=m(q+i);               %定义一个1*p的向量
end
B=[-r',R];                                  %对应矩阵B
[U,K,V]=svd(B);              %求增广矩阵B的SVD，并存储奇异值和矩阵V
P=rank(B);                  % 求得增广矩阵B的有效秩，定义为P
S=zeros(P+1);          %构造一个(p+1)*(p+1)维的矩阵S                   

for j=1:p                  
for i=1:p+1-P
  djj=K(j,j)*K(j,j);     %构造djj,并求其平方
  vij=V(i:i+p,j);         %构造vij
 S= S+djj*vij*vij';     %计算矩阵S的二重级数求和
 end
end

Sni=inv(S);                   %求S逆矩阵
a=zeros(1,P);                     %构造a矩阵
for i=1:P                                
 a(1,i)=Sni(i+1,1)/Sni(1,1);     %求出矩阵AR模型参数a=[a1,...,ap]'
end

a1=fliplr(a)        %将a进行元素对调,使a1=[ap,...,a1]'
                                         

figure(1); 
freqz(1,a1,128,1);            %求出信号的幅频响应和相频响应波形
title('AR模型的总体最小二乘（SVD-TLS）估计');
legend(strcat('AR（p）阶数p=',int2str(p)));
grid on;

这时我们取的矩阵A阶数为4，当然可以取其他更大的值，来得到最后的图像：

这里大致介绍了SVD-TLS的具体实现步骤，但是由于第一次在网上总结，这个编辑器也不大会用，所以准备了很多推导内容也没有呈现出来。