现代功率谱估计(3):SVD-TLS,奇异值分解—总体最小二乘方法求解AR模型参数
Yuler-Walker方程及修正Yuler-Walker方程
对于一个AR过程,其输出信号的自相关函数和AR系数有以下方程:

现在的问题是,如何求解AR模型的系数以及AR模型的阶数,根据AR模型的参数可辨识性定理:

可以通过求解有限个(个)修正Yuler-Walker方程来进行求解AR模型的系数,接下来的问题是:如何确定AR模型的合适的阶数,从而得知上面方程的个数?
为了解决这一问题我们先假设一个较大的AR模型的阶数,在假设阶数为的情况下,上面的修正Yuler-Walker方程依旧成立:

其中上面在我们假设,的情况下,需要满足以下命题:

接下来就可以使用奇异值分解的方法,利用上面矩阵的奇异值来确定AR模型的真实阶数
到这里先停一下,先放上具体的代码示例,动手将上面矩阵构建出来:
%% 根据样本序列估计的自相关生成矩阵Re
Rxx=xcorr(xn);
pe=100; % pe>=p;
qe=150; % qe>=q,(qe-pe)>=(q-p);
M=300; % M>=pe;
for i = 1 : N -1
Rx(i) = Rxx(N - 1 + i); % xcorr函数默认的自相关函数为最大值的点并非位于横坐标0处,此处将其平移到坐标轴处
end
Re=zeros(M,pe+1); % 书中p129, 4.4.23
for i=1:M
for j=1:(pe+1)
Re(i,j)=Rx(qe + i - j + 1);
end
end
上面生成的矩阵Re即为根据样本估计自相关得到的修正Yuler-Walker方程系数矩阵。
奇异值分解确定AR模型的阶数
先对上面得到的样本估计相关函数矩阵做奇异值分解
再使用归一化奇异值的比值,设定合适的阈值来确定上面矩阵的有效秩,归一化奇异值定义为:

通过选择合适的阈值(一个接近于0的正数),并将大于此阈值的最大整数(奇异值个数)取作矩阵的有效秩,即可得到一个相比于之前设定的更加合适的AR模型。
这一段通过奇异值分解,根据奇异值的大小比较确定AR模型的阶数的代码如下:
%% 根据前几个sigma的归一化比值,确定AR模型的阶数p
[U,S,V]=svd(Re); % 奇异值分解
%%求p值
p=0;
SS=zeros(1,pe+1)';
for i=1:(pe+1)
if S(i,i)/S(1,1)>0.05 % p132页4.4.40, 归一化比值,挑出前几个较大的奇异值,根据其个数确定阶数p
p=p+1;
if p>0
SS(p)=S(i,i);
end
end
end
总体最小二乘法(TLS)确定AR模型的参数
上面我们通过随机信号的样本信号估计出了修正Yuler-Walker方程的系数矩阵,通过奇异值分解,归一化比值和阈值比较之后确定了AR模型的阶数,接下来的问题就是如何求解修正Yuler-Walker方程中的未知数:

根据张贤达现代信号处理书上p135,式4.4.56,式4.4.58步骤:

其中为矩阵第个奇异值的平方,向量为矩阵奇异值分解之后的酉矩阵的第列的一个加窗段:

构造出矩阵之后,即可求出上面解修正Yuler-Walker方程中的未知数。具体计算公式为:

这一段根据矩阵奇异值分解之后的酉矩阵和相应奇异值构建矩阵的代码如下:
[U, S, V] = svd( Re ); % 奇异值分解
%%求p值
p = 0;
SS=zeros(1, pe + 1)';
for i = 1 : ( pe + 1 )
if S(i , i) / S(1 , 1) > 0.05 % p132页4.4.40, 归一化比值,挑出前几个较大的奇异值,根据其个数确定阶数p
p = p + 1;
if p > 0
SS( p ) = S(i , i );
end
end
end
%% 求Sp及其逆
Sp = zeros( p + 1 , p + 1 );
for j = 1 : p
for i = 1 : (pe+1-p)
Sp = Sp + SS(j)^2 * V(i:i+p , j) * V( i : i+p , j )'; % p135 式4.4.58
end
end
至此,AR模型阶数已确定,参数也已求出,接下来可利用Cadzow谱估计子对功率谱进行估计(这部分代码以后再补上)
最后,给出汇总的实验仿真代码:
clc;
clear;
close all;
%% 参数设置
f1 = 70; f2 = 100; f3 = 120;
N = 1000;
fs = 512;
n = (1: 1000) * 1/ fs;
%% 仿真数据产生
xn = 2 * sin(2 * pi * f1 * n + pi / 3) + 2 * sin(2 * pi * f2 * n + pi / 4) + 2 * sin(2 * pi * f3 * n + pi / 5) ;
xn = xn + 2 * randn(size(n));
%% 根据样本序列估计的自相关生成矩阵Re
Rxx = xcorr(xn);
pe = 100; % pe>=p;
qe = 150; % qe>=q,(qe-pe)>=(q-p);
M = 300; % M>=pe;
for i = 1 : N -1
Rx(i) = Rxx(N - 1 + i); % xcorr函数默认的自相关函数为最大值的点并非位于横坐标0处,此处将其平移到坐标轴处
end
Re=zeros(M, pe+1); % 书中p129, 4.4.23
for i=1: M
for j=1: (pe+1)
Re(i, j)=Rx(qe + i - j + 1);
end
end
%% 根据前几个sigma的归一化比值,确定AR模型的阶数p
[U, S, V] = svd( Re ); % 奇异值分解
%%求p值
p = 0;
SS=zeros(1, pe + 1)';
for i = 1 : ( pe + 1 )
if S(i , i) / S(1 , 1) > 0.05 % p132页4.4.40, 归一化比值,挑出前几个较大的奇异值,根据其个数确定阶数p
p = p + 1;
if p > 0
SS( p ) = S(i , i );
end
end
end
%% 求Sp及其逆
Sp = zeros( p + 1 , p + 1 );
for j = 1 : p
for i = 1 : (pe+1-p)
Sp = Sp + SS(j)^2 * V(i:i+p , j) * V( i : i+p , j )'; % p135 式4.4.58
end
end
Sp_=inv(Sp);
%% 总体最小二乘估计
x=zeros(1,p+1)';
for i=1:(p+1)
x(i)=Sp_(i, 1)/ Sp_(1, 1); % p135 4.4.58
end
%% 求响应
a=x;
%% 随机信号的幅度谱
subplot 221
xfft = fft(xn);
freq = (1: length(xfft))/length(xfft) * fs;
f = freq(1:(length(xfft)/2 + 1));
plot(f, abs(xfft(1 : (length(xfft)/2)+1)));
title("原始信号的幅度谱")
%% AR模型的功率谱密度分布情况
[H, w] = freqz(1, a, fs);
subplot 222
plot(w/(2 * pi) * fs, abs(H).^2);
title("p阶AR模型的功率谱密度分布情况")
%% 周期图法
subplot 223
window=boxcar(length(xn));
nfft=fs;
[Pxx, w]=periodogram(xn,window,nfft,fs); %直接法
plot(w, Pxx)
title("直接周期图法计算随机信号Xn的功率谱密度分布")
运行结果:

【推荐】编程新体验,更懂你的AI,立即体验豆包MarsCode编程助手
【推荐】凌霞软件回馈社区,博客园 & 1Panel & Halo 联合会员上线
【推荐】抖音旗下AI助手豆包,你的智能百科全书,全免费不限次数
【推荐】博客园社区专享云产品让利特惠,阿里云新客6.5折上折
【推荐】轻量又高性能的 SSH 工具 IShell:AI 加持,快人一步
· 为DeepSeek添加本地知识库
· 精选4款基于.NET开源、功能强大的通讯调试工具
· DeepSeek智能编程
· 大模型工具KTransformer的安装
· [计算机/硬件/GPU] 显卡