现代谱估计实验报告

杭州电子科技大学现代谱估计实验报告

现代谱估计实验报告

吴迪松 20040089 信息与信息处理

一、实验题目：

用MATLAB 编写一下算法的程序：

(1)修正协方差法

(2)多重信号分类（MUSIC ）算法

(3)ESPRIT 算法

(4)皮萨论科（Pisarenko ）谐波分解法

并对各算法进行分析。

二、实验原理

(1) 修正协方差法

1、得出输入信号x(n),此输入信号是与白噪声相加的。

xx R ˆ()()()x

n x n n t =+ 2、求自相关估计器(),xx c j k

()()11**01,()(2()N N P xx n p n c j k x n j x n k x n j x n k N p −−−==)()⎡⎤=−−+++⎢⎥−⎣⎦

∑∑ 3、求滤波器的系数 ˆ()a

k

1ˆ()(,)(,1)p xx xx k a

k c l k c l ==−∑ 4、求白噪声方差的估计值2ˆσ

2

1

ˆˆ(1,1)()(1,)p xx xx k c a k c σ==+∑k 5、求功率谱

()p f 2221ˆ()ˆ1()p j fk k p f a

k e πσ

−==+∑

(2)多重信号分类（MUSIC ）算法

1、估计样本自相关矩阵 ˆxx r

2、对作特征值分解 ˆxx r

3、确定的最小特征值分解的数目,求出这个最 ˆxx r

E n E n 小特征值1,,p M λλ+"，令

()2121p p M E

n σλλλ++=+++" 并求出相对应的特征向量,构造噪声矩阵

1,,p v v +"M

21,,p M V v v +⎡⎤=⎣⎦"

4、构造函数：

()211

ˆMU M H i i p P f e v =+=∑

(3)ESPRIT 算法

1、利用已知信息求 {}(0),(1),()xx xx xx r r r M "

2、由{}(0),(1),()xx xx xx r r r M "构造自相关矩阵xx R 和互相关矩

阵xy R

3、计算xx R 的特征值分解。对于M>P ，最小特征值为

噪声方差2σ

4、利用2σ计算2xx xx C R I σ=−和2xy xy C R I σ=−

5、计算矩阵对{},xx xy C C 的广义特征值分解。从而确定谐

波频率。

(4)皮萨论科（Pisarenko ）谐波分解法

1、求x(n)的自相关函数

xx r 2、求的toeplitz 矩阵xx r xx R

3、求出xx R 的特征值，由此得出最小特征值

4、求出最小特征值对应的一列特征向量

5、求这个特征矢量形成的多项式的根。由这些根求出

它的角度，并且最后求出频率

三、各算法的MATLAB 程序

%AR model parameter estimate

clear

m=sqrt(-1);

delta=0.1;

a1=0.5;

sample=32; %number of sample spot

p=10; %number of sample spot in coef method;

f1=0.05; f2=0.20; f3=0.45; %三个频率分量

fstep=0.01;

fstart=-0.5;

fend=0.5;

f=fstart:fstep:fend; %频率的范围

nfft=(fend-fstart)/fstep+1;

%un=urn+juin

urn= normrnd(0,delta/2,1,sample);

uin= normrnd(0,delta/2,1,sample);

un=urn+m*uin;

%calculate zn 噪声

for n=1:sample-1

zn(1)=un(1);

zn(n+1)=a1*zn(n)+un(n+1);

end

%calculate xn 出入信号

for n=1:sample

xn(n)=2*cos(2*pi*f1*(n-1))+2*cos(2*pi*f2*(n-1))+2*cos(2*pi*f3*(n-1)) +sqrt(2)*real(zn(n));

end %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%calculate cxx(j,k) 公式3.102

for j=1:p

for k=1:p

s=0;

for n=p+1:sample

s=s+1/2/(sample-p)*(conj(xn(n-j))*xn(n-k)+xn(n-p+j)*conj(xn(n-p+k)));

end

cxx(j,k)=s;

end

end

%calculate cxx0(j,1) 公式 3.103

for j=1:p

s=0;

for n=p+1:sample

s=s+1/2/(sample-p)*(conj(xn(n-j))*xn(n)+xn(n-p+j)*conj(xn(n-p))); end

cxx0(j,1)=s;

end

%calculate a 公式3.104

a=-inv(cxx)*cxx0;

%a(k)*exp(-j2pifk)累加

for i=1:length(f)

sum=0;

for k=1:p

sum=sum+a(k)*exp(-m*2*pi*f(i)*k);

end