功率谱估计Levinson递推法和Burg法..docx
《功率谱估计Levinson递推法和Burg法..docx》由会员分享，可在线阅读，更多相关《功率谱估计Levinson递推法和Burg法..docx（37页珍藏版）》请在冰点文库上搜索。
数字信号处理实验报告
姓名：
学号：
日期：
2015.12.14
1.实验任务
信号为两个正弦信号加高斯白噪声，各正弦信号的信噪比均为10dB，长度为，信号频率分别为和，初始相位，取，取不同的数值：
0.3，0.25。
为采样率。
（1）分别用Levinson递推法和Burg法进行功率谱估计，并分析改变数据长度、模型阶数对谱估计结果的影响。
（2）当正弦信号相位、频率、信噪比改变后，上述谱估计的结果有何变化？
并作分析说明。
2.原理分析
2.1现代谱估计中的参数建模
根据参数模型来描述随机信号的方法，我们可以知道，如果能确定信号的信号模型，根据信号观测数据求出模型参数，系统函数用表示，模型输入白噪声，其方差为，信号的功率谱用下式求出：
按照这种求功率谱的思路，功率谱估计可分为三个步骤：
（1）选择合适的信号模型；
（2）根据有限的观测数据，或者它的有限个自相关函数的估计值，估计模型的参数；
（3）计算墨香的输出功率谱。
其中以
（1）、
（2）两步最为关键。
按照模型的不同，谱估计的方法有许多种，它们共同的特点是对信号观测区以外的数据不假设为0，而先根据信号观测数据估计模型参数，按照求模型输出功率的方法估计信号功率谱，回避了数据观测区以外的数据假设问题。
下面分析AR谱估计的两种方法：
自相关法——列文森（Levenson）递推法和伯格（Burg）递推法。
这两种方法均为已知信号观测数据，估计功率谱，两者共同特点是由信号观测数据求模型系数时采用信号预测误差最小的原则。
对于长记录数据，这些方法的估计质量是相似的，但对于短记录数据，不同方法之间存在差别。
2.2自相关法——列文森（Levenson）递推法
自相关法的出发点是选择AR模型参数使预测误差功率最小，预测误差功率为
假设信号的数据区在范围，有个预测系数，个数据经过冲激响应为的滤波器，输出预测误差的长度为，因此应用下式计算：
的长度长于数据的长度，上式中数据的两端需补充零点，相当于对无穷长的信号加窗处理，得到长度为N的数据。
上式对系数的实部和虚部求微分使预测误差功率最小，得到
（1）
式中自相关函数采用有偏自相关估计，即
对比上式，可知式
（1）即为已推导出的Yule-Walker方程，因此自相关法也是基于解Yule-Walker方程的一种方法。
但是直接解该方程，需要计算逆矩阵，不方便，因此，基于Yule-Walker方程中自相关矩阵的性质，导出Levinson-Durbin递推法，这是一种高效的解方程的方法。
Levinson-Durbin算法首先由一阶AR模型开始：
一阶AR模型的Yule-Walker方程为
由该方程解出
然后令，以此类推，可以得到一般递推公式如下：
称为反射系数，。
，随着阶数增加，预测误差功率将减少或不变。
由k=1开始递推，递推到k=p，依次得到各阶模型参数，
AR模型的各个系数及模型输入白噪声方差求出后，信号功率谱用下式计算
这种方法计算简单，但需要预先估计出信号自相关函数，实际中只能按照信号的有限个观测数据估计自相关函数。
当观测数据长度较短时，估计误差较大，会出现谱峰频率偏移和谱线分裂（在信号谱峰附近产生虚假谱线）；如数据很长，估计自相关函数较准确，但计算量大，应适当选择数据长度。
2.3伯格（Burg）递推法
Levinson-Durbin递推法需要由观测数据估计自相关函数，这是它的缺点。
而伯格递推法则由信号观测数据直接计算AR模型参数。
伯格递推法利用Levinson-Durbin递推公式，导出前向预测误差与后向预测误差，并按照使它们最小的原则求出，从而实现不用估计自相关函数，直接用观测数据得出结果。
Burg递推法思想：
借助格型预测误差滤波器，求前向、后向预测误差平均功率，选择使其最小，求出。
之后，再利用Levinson-Durbin递推法求模型参数和输入噪声方差。
设信号的观测数据区间：
前向、后向预测误差功率分别用和表示，预测误差平均功率用表示，公式分别为
前向、后向观测误差公式分别为
上式中，信号项的自变量最大的是n,最小的是n-p，为了保证计算范围不超出给定的数据范围，在和计算公式中，选择求和范围为：
。
为求预测误差平均功率最小时的反射系数，令,将前、后向预测误差的递推公式代入得
Burg递推法求AR模型参数的递推公式总结如下：
（1）
（2）
（3）
（4）
（5）
（6）
（7）
3.编程思想
（1）编写程序产生题目要求的信号和噪声
（2）然后分别用两种方法的递推流程进行谱估计
（3）改变题目中要求的变量参数，分析结果的变化
4.代码
Levenson
clc;
clearall;
fs=100;%采样频率
Ts=1/fs;
N=2^7;%数据长度
p1=20;%阶数
f1=0.2*fs;f2=0.25*fs;%设置信号频率
pha1=0;pha2=0;%初始相位
SNR=2;%设置信噪比
%产生信号
w=randn（1,N）;
Am=sqrt（2*10^（SNR/10））;
k=1:
N;
x1=Am*sin（2*pi*f1*k*Ts+pha1）;x2=Am*sin（2*pi*f2*k*Ts+pha2）;
xx=x1+x2;x=w+x1+x2;
figure
（1）
subplot（2,1,1）;plot（xx）;ylim（[-20,20]）;title（'两个正弦信号波形'）;gridon;
subplot（2,1,2）;plot（x）;ylim（[-20,20]）;title（'加噪信号波形'）;gridon;
%计算自相关函数
R=[];
form=1:
N
s=0;
forn=1:
N-（m-1）
s=s+x（n）*x（n+m-1）;
end
R（m）=s/N;
end
%列文森递推法
a（1,1）=-R
（2）/R
（1）;cov
（1）=R
（1）+a（1,1）*R
（2）;
forp=2:
p1
sum=0;
fork=1:
p-1
sum=sum+a（p-1,k）*R（p-k+1）;
end
a（p,p）=-（R（p+1）+sum）/cov（p-1）;%计算反射系数
fork=1:
p-1
a（p,k）=a（p-1,k）+a（p,p）*a（p-1,p-k）;
end
cov（p）=（1-（abs（a（k,k）））^2）*cov（p-1）;
end
%计算功率谱，功率谱以2*pi为周期，又信号为实信号，只需输出0到P1即可；
W=0.01:
0.01:
pi;Z=0;
fork=1:
p1
Z=Z+（a（p1,k）.*exp（-j*k*W））;
end
PSD=1./（（abs（1+Z））.^2）;%算出功率谱
F=W*fs/（pi*2）;%角频率坐标换算成频率坐标
figure
（2）
plot（F,abs（PSD））;xlabel（'频率（Hz）'）;ylabel（'功率谱'）;
title（'自相关法-列文森Levenson递推法的功率谱估计'）;grid;
figure（3）
p=1:
p1;plot（p,cov（p）,'b'）;xlabel（'模型阶数'）;ylabel（'预测误差功率大小'）;
title（'预测误差功率变化趋势'）;grid;
Burg
clc;
clearall;
fs=100;%采样频率
Ts=1/fs;
N=2^8;%数据长度
p1=20;%阶数
f1=0.2*fs;
f2=0.25*fs;%设置信号频率
pha1=0;pha2=0;
SNR=2;%设置信噪比
%产生信号
w=randn（1,N）;%噪声为高斯白噪声，功率为1.
Am=sqrt（2*10^（SNR/10））;
k=1:
N;
x1=Am*sin（2*pi*f1*k*Ts+pha1）;x2=Am*sin（2*pi*f2*k*Ts+pha2）;
x=w+x1+x2;
%递推
ef=zeros（p1,N）;eb=zeros（p1,N）;
ef（1,:
）=x;
eb（1,:
）=x;
cov
（1）=x*x'/N;
k
（1）=0;
a=zeros（p1+1,p1+1）;
forp=2:
p1+1
numerator=0;
denominator=0;
forn=p:
N
numerator=numerator+（-2）*ef（p-1,n）*eb（p-1,n-1）;
denominator=denominator+（ef（p-1,n））^2+（eb（p-1,n-1））^2;
end
k（p）=numerator./（denominator+0.0001）;
cov（p）=（1-abs（k（p））^2）.*cov（p-1）;
a（p,p）=k（p）;
forn=p:
N
ef（p,n）=ef（p-1,n）+k（p）.*eb（p-1,n-1）;
eb（p,n）=eb（p-1,n-1）+k（p）.*ef（p-1,n）;
end
end
k=k（1,2:
p1+1）;a=a（2:
p1+1,2:
p1+1）;
forp=2:
p1
fori=1:
p-1
a（p,i）=a（p-1,i）+k（p）*a（p-1,p-i）;
end
end
%功率谱
Z=0;W=1:
0.01:
pi;
fori=1:
p
Z=Z+（a（p,i）*exp（-j*W*i））;
end;
pxx=cov（p1+1）./（abs（1+Z）.^2）;
F=W*fs/（pi*2）;%角频率坐标换算成频率坐标
figure
（1）
plot（F,pxx）;
xlabel（'Frequency（Hz）'）;
ylabel（'PowerSpectralDensity'）;
title（'伯格（Burg）递推法的功率谱估计'）;
grid;
%figure
（2）
%p=1:
p1;
%plot（p,cov（p）,'b'）;
%xlabel（'模型阶次'）;
%ylabel（'预测误差功率大小'）;
%title（'预测误差功率的变化趋势'）;
%grid;
5.实验结果及分析
5.1产生信号
5.2Levenson递推法
1.数据长度的影响（阶数设置为20阶）
N=32
N=64
N=128
N=1024
N=8192
N=16384
2.模型阶数的影响（数据长度设置为N=1024）
预测误差功率的变化（从1阶到50阶）
不同阶数（p1表示）功率谱的图像如下：
P1=5
P1=10
P1=20
P1=30
P1=50
3.相位的影响（设置N=128p1=20）
0
Pi/6
Pi/3
Pi/2
Pi
4.频率的影响（N=128p1=20初始相位为0）
Fs=100,f1=0.2*fsf2=0.25*fs
Fs=100,f1=0.2*fsf2=0.3*fs
Fs=1000,f1=0.2*fsf2=0.25*fs
Fs=1000,f1=0.2*fsf2=0.3*fs
5.信噪比的影响
SNR=20dB
SNR=10dB
SNR=6dB
SNR=2dB
5.3Burg递推法
1.数据长度的影响
N=32
N=64
N=256
N=1024
N=8192
N=16384
2.模型阶数的影响（N=128）
预测误差功率的变化趋势
不同阶数的功率谱
P1=5
P1=10
P1=20
P1=50
3.相位的影响（N=256p1=20）
0
Pi/6
Pi/3
Pi/2
Pi
4.频率的影响
Fs=100,f1=0.2*fsf2=0.25*fs
Fs=100,f1=0.2*fsf2=0.3*fs
Fs=1000,f1=0.2*fsf2=0.25*fs
Fs=1000,f1=0.2*fsf2=0.3*fs
分析：
5.信噪比的影响
SNR=20dB
SNR=10dB
SNR=6dB
SNR=2dB
分析：
信噪比越大，功率谱估计效果越好。
信噪比较小的时候，可以看到，谱估计图像中，频率在20Hz和25Hz的地方，频谱不够尖锐，
升级会员 