关于SVD

萌宠菠菠乐园
2024-10-04 01:46

SVD-TLS算法

在许多工程应用中，功率谱的分析和估计是很重要的，因为它能够给出被分析对象的能量随频率的变化情况。在识别一个目标时，功率谱可以当做一个目标特征。
估计功率谱密度的平滑周期图是一类非参数化方法，与任何模型参数无关。它的主要问题是：由于假定信号的自相关函数在数据观测区以外为零，因此估计出来的功率谱很难与信号的真实功率谱相匹配。
而SVD-TLS是一个在现代谱估计中一个基础的算法。
SVD代表的是奇异值分解，它主要用于求解线性方程组，对于奇异值的使用和证明在现代信号处理的现代谱中描述很详细，最后可表明矩阵A(K)逼近A的精确度取决于被置零的那些奇异值的平方和。
对于k取到某一个合适值p时，逼近误差向量A-A(k)的Frobenious范数已足够小，并且当k>p时，该范数不再明显减小。这一p值称为矩阵A的有效秩。确定矩阵A的有效秩的这一方法可用以下归一化比值法：
归一化比值定义
这时需要先确定一个接近1的阈值，用来确定矩阵A的有效秩。
同时还有一个方法就是归一化奇异值来确定有效秩：
归一化奇异值法
与上面的归一化比值相反，需要确定一个接近零的正整数作为阈值来确定矩阵A的有效秩p.
TLS也就是最小二乘法，通过奇异值分解我们可以得到一个矩阵A的有效秩p，这时可以利用最小二乘法确定p个AR参数的估计值。这里有两个问题：

必须重新列出法方程组，使它只包含p个未知数；求解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

12345678910111213141516171819

我们也可以带入时间序列确定其自相关函数，构造一个矩阵来运用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;

12345678910111213141516171819202122232425262728293031323334353637383940414243444546474849505152535455565758

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