现代谱估计
现代谱估计实验报告
1 实验目的
功率谱估计在实际工程中有重要应用价值。如在语音信号识别、雷达杂波分析、波达方向估计、地震勘探信号处理、水声信号处理、系统辨识中非线性系统识别、物理光学中透镜干涉、流体力学的内波分析、太阳黑子活动周期研究等许多领域发挥了重要作用。
本次实验的目的主要是深入理解现代谱估计的基本理论，包括ARMA 模型、ARMA 谱估计。掌握现代谱估计的基本方法，包括SVD-TLS 算法等。利用ARMA 功率谱估计中Cadzow 谱估计子和Kaveh 谱估计子来进行谱估计。
2 实验原理
2.1 背景
若离散随机过程{x(n)}服从线性差分方程
x (n ) +∑a i x (n -i ) =e (n ) +∑b j e (n -j )
i =1j =1p q （1）
式中e （n ）是一离散白噪声，则称{x(n)}为ARMA 过程，而式（1）所示的差分方程称为ARMA 模型。系数a 1,a 2……a p , 和b 1,b 2……b q , 分别称为自回归参数和滑动平均参数，而p 和q 分别叫做AR 阶数和MA 阶数。式（1）所示的ARMA 过程，其功率谱密度为
2P x (w ) =δB (z ) A (z ) z ==δ2e jw B (e ) B (e ) jw jw （2）
ARMA 谱估计的目的是使用N 个已知的观测数据x(0)，x(1)…..x(N-1)计算出ARMA 过程{x(n)}的功率谱密度估计。
在实际中，可以运用cadzow 谱估计子和kaveh 谱估计子来估计，cadzow 谱估计子秩序确定AR 阶数p 和估计AR 参数，而kaveh 谱估计子也只需要确定AR 阶数p 和估计AR 参数以及MA 阶数。
2.2 相关算法
AR 阶数p 的确定用奇异值分解（SVD ）,AR 参数的估计用总体最小二乘法(TLS),即应用（SVD —TLS ）算法来完成ARMA 谱估计。
SVD —TLS 算法：
步骤1 计算增广矩阵B 的SVD ，并储存奇异值和矩阵V;
步骤2 确定增广矩阵B 的有效秩p ；
步骤3 计算矩阵S;
步骤4 求S 的逆矩阵S--, 并计算出未知参数的总体最小二乘估计。 3 实验内容
仿真的观测数据由下式给出：
（3） xn = square(W*n)+0.2*randn(1,N)
其中，fs = 20000，n = 0:1/fs:0.1，N = length(n)，W = 2000*pi。
1、采样周期图法进行谱估计
2、假设AR 阶数未知，用SVD-TLS 方法确定AR 阶数和参数，然后使用Cadzow 谱估计子进行谱估计。
4 Matlab仿真
仿真的观测数据时域信号如图1所示。
图1 观测数据时域信号
1、经典功率谱估计
周期图法是把随机序列x(n)的N 个观测数据视为一能量有限的序列。直接计算x(n)的离散傅立叶变换，得X(k)，然后再取其幅值的平方，并除以N ，作为序列x(n)真实功率谱的估计。仿真图如图2所示。
图2 周期图法功率谱
2、现代功率谱估计
现代功率谱估计即参数谱估计方法是通过观测数据估计参数模型再按照求参数模型输出功率的方法估计信号功率谱。主要是针对经典谱估计的分辨率低和方差性能不好等问题提出的。按照上面介绍的步骤，编写程序对观测信号x(n)进行仿真，可以设置不同的M ，qe ，pe 的值，以便分析对比。图3是设置了M=2001，qe=100，pe=50，后得出的x(n)的功率谱图形。
图3 ARMA模型功率谱
5 实验总结
本次实验分别用了周期图法和ARMA 模型的参数估计方法对方波信号进行了功率谱估计，通过实验和得到的仿真图对比可以发现：
通过周期图法得到的功率谱估计频谱分辨率较低，不能适应高分辨率功率谱估计的要求，参数化的谱估计可以获得高频率分辨率的功率谱。经典功率谱估计的分辨率反比于有效信号的长度，但现代谱估计的分辨率可以不受此限制。这是因为对于给定的N 点有限长序列x(n)，虽然其估计出的自相关函数也是有限长的， 但是现代谱估计的一些隐含着数据和自相关函数的外推，使其可能的长度超过给定的长度，不像经典谱估计那样受窗函数的影响。因而现代谱的分别率比较高，而且现代谱线要平滑得多，从上图可以清楚看出。
6 附录
Matlab 程序如下：
main.m
clear;
close all
fs = 20000;
n = 0:1/fs:0.1;
N = length(n);
W = 2000*pi;
x1n = square(W*n);
x2n = randn(1,N);
xn = x1n+0.2*x2n;
figure;plot(n,xn);
title(' 时域信号' );
Nfft = 100;
[Pxx,f] = period(xn,fs,Nfft);
figure;plot(f,Pxx);
title(' 周期图法功率谱' );
%ARMA谱估计
pe = 50;
qe = 100;
NARMA = length(xn);
M = length(n);
[a,Rx,p] = ARMA (xn,qe,pe,M);
%Cadzow谱估计子
[Pw] = Cadzow(a,Rx,p,NARMA);
%功率谱×
figure;plot((0:length(Pw)-1)*fs/length(Pw),Pw);
title('ARMA 模型' );
period.m
function [Pxx,f] = period(xn,fs,Nfft)
Pxx = abs(fft(xn,Nfft).^2)/Nfft;
f = (0:length(Pxx)-1)*fs/length(Pxx);
ARMA.m
function [a,Rxx,p] = ARMA(xn,qe,pe,M)
Rxx = xcorr(xn,'unbiased' );
for (i = 1:M)
for (j = 1:pe+1)
Re(i,j) = Rxx(pe+i+1-j);
end
end
[U,S,V] = svd(Re);
Ak = 0;
for (i = 1:pe+1)
Ak = Ak + S(i,i)^2;
end ;
Akf = 0;
v = Akf/Ak;
p = 0;
while v
p = p+1;
Akf = Akf + S(p,p)^2;
v = Akf/Ak;
end
Sp = 0;
for j =1:p
for i =1:pe+1-p
Sp=Sp+S(i,i)^2*V(i:i+p,j)*V(i:i+p,j)';
end
end
invSp = inv(Sp);
for (i = 1:p)
a(i) = invSp(i+1,1)/invSp(1,1);
end
Cadzow.m
function [Pw] = Cadzow(a,Rx,p,N)
Ro = Rx;
Ro(N) = 0.5*Rx(N);
p = p-1;
n = zeros(1,p+1);
for k=1:p+1
n(k) = a(1:p+1)*Ro(k+N:-1:k-p+N)';
end
wfreq = linspace(0,2*pi,N);
z = exp(sqrt(-1)*wfreq);
Pw = polyval(n,z)./polyval(a,z)+polyval(n,z.^-1)./polyval(a,z.^-1); Pw = abs(Pw);