现代谱估计

现代谱估计实验报告现代谱估计实验报告 1 1 实验目的实验目的 功率谱估计在实际工程中有重要应用价值。如在语音信号识别、雷达杂波分 析、波达方向估计、地震勘探信号处理、水声信号处理、系统辨识中非线性系统 识别、物理光学中透镜干涉、流体力学的内波分析、太阳黑子活动周期研究等许 多领域发挥了重要作用。 本次实验的目的主要是深入理解现代谱估计的基本理论，包括 ARMA 模型、 ARMA 谱估计。 掌握现代谱估计的基本方法， 包括 SVD-TLS 算法等。 利用 ARMA 功率谱估计中 Cadzow 谱估计子和 Kaveh 谱估计子来进行谱估计。 2 2 实验原理实验原理 2.1 背景 若离散随机过程{x(n)}服从线性差分方程 x(n)   a i x(n  i)  e(n)   b je(n  j) i1j 1 pq （1） 式中 e（n）是一离散白噪声，则称{x(n)}为 ARMA 过程，而式（1）所示的差分 方程称为 ARMA 模型。 系数 a1,a2……ap,和 b1,b2……bq,分别称为自回归参数和滑 动平均参数，而p 和 q 分别叫做 AR 阶数和 MA 阶数。式（1）所示的ARMA 过 程，其功率谱密度为 2P x (w)   B(z) A(z) z  2 ejw 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 实验内容实验内容 仿真的观测数据由下式给出： （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 Matlab4 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 5 实验总结实验总结 本次实验分别用了周期图法和 ARMA 模型的参数估计方法对方波信号进行 了功率谱估计，通过实验和得到的仿真图对比可以发现： 通过周期图法得到的功率谱估计频谱分辨率较低，不能适应高分辨率功率谱 估计的要求， 参数化的谱估计可以获得高频率分辨率的功率谱。经典功率谱估计 的分辨率反比于有效信号的长度，但现代谱估计的分辨率可以不受此限制。这是 因为对于给定的 N 点有限长序列 x(n)，虽然其估计出的自相关函数也是有限长 的， 但是现代谱估计的一些隐含着数据和自相关函数的外推，使其可能的长度 超过给定的长度， 不像经典谱估计那样受窗函数的影响。因而现代谱的分别率比 较高，而且现代谱线要平滑得多，从上图可以清楚看出。 6 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 0.99995 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(