谱估计实验报告汇总文档格式.docx
《谱估计实验报告汇总文档格式.docx》由会员分享，可在线阅读，更多相关《谱估计实验报告汇总文档格式.docx（23页珍藏版）》请在冰点文库上搜索。
fori=2:
N
x（i）=0.8*x（i-1）+w（i）;
end
此方法简单，可以很容易地产生所需数目的数据。
第二种思路：
利用卷积的方法
对线性时不变系统，输入输出满足卷积关系：
，可得
，从而可得系统的冲击响应：
然后进行卷积运算即可。
Matlab代码实现如下：
n=1:
N;
h（n）=（0.8）.^（n-1）;
%要注意n-1不是n，因为冲击响应是从0开始的
y=conv（w（n）,h（n））；
%共2*N-1项，取前N项即可
需注意：
实际h（n）是从0开始的，matlab处理元素从下标1开始，因此，公式中应是n-1不是n。
而且，计算完成后卷积结果是为2*N-1项，取前N项即可。
两种方法结果
为方便观察，令N=5时，实验结果如下：
x=
0.62321.29761.97900.59110.6849
y=
0.62321.29761.97900.59110.68490.34370.0131-0.29780.0868
取卷积的前N项，可以看出两种方式结果是相同的。
（2）、信号x（n）的理想功率谱
系统为AR模型，理想功率谱为：
因此，对h（n）进行傅里叶变换后，取模的平方即可。
（3）、周期图法谱估计
根据周期图法谱估计原理：
先对观测数据x（n）进行傅里叶变换，后平方，最后除以N即可。
（4）、Bartlett平均周期图谱估计
为了减小估计的方差，将数据分为L段，则每段有M=N/L个数据，分别用周期图法进行估计后求平均。
具体公式如下：
将得到的L个周期图进行平均，作为信号x（n）的功率谱估计，公式如下：
经过平均处理，可使估计方差减小。
2、实验结果与分析
图1白噪声w（n）
图2平稳随机信号x（n）
信号理想功率谱如下图3所示：
图3信号的理想功率谱
从图中可以看出，理想功率谱是平滑的。
下图4是功率谱的分贝形式：
图4信号的理想功率谱（dB）
当数据长度N为1024时，实验结果如下图5所示：
图5N=1024周期图谱估计结果
当数据长度N为256时，实验结果如下图6所示：
图6N=256周期图谱估计结果
对比图如下图7：
图7N=1024/256周期图谱估计对比
从以上结果可以看出N=1024时频谱分辨率明显要高于N=256时的频谱分辨率。
当数据长度N=1024及256时，分段数L分别为2和8时信号
的功率谱为便于对比，将结果表示如下图8一幅图中：
图8Bartlett平均周期图谱估计
结果分析：
1、首先数据的长度对分辨率有影响，数据长度N=1024时的频谱分辨率比N=256时的频谱分辨率高；
2、分段数L对频谱的分辨率和平滑性（方差）也有很大影响。
当数据数目N一定时，L加大，每一段的数据量M就会减小，因此估计方差减小，曲线越平滑，但此时偏移加大，分辨率降低，即估计量的方差和分辨率是一对矛盾，二者的效果可以通过合理选取L互换。
3、收获体会
1、通过实验，对经典谱估计法有了更深刻的理解，根据课堂中经典谱估计的理论，掌握了周期图法，分段周期图法的具体实现，更加深刻地体会到了理论的原理以及这些估计方法的优缺点，对这些估计方法获得了真正的理解；
2、很小的细节也要注意，也需要认真思考，如信号x（n）的产生，通过两种方法的对比，对卷积有了更深入的认识，通过不同方法得到相同结果，学会了从不同的角度分析问题；
3、对matlab掌握更好，如自己利用子函数调用的方式实现了对不同分段平均周期图法的实现。
函数传递参数分段数L即数据x（n），子函数代码如下：
functionPx_L=Px_mean（L,x）
N=length（x）;
%计算数据长度
M=N/L;
%每段有M个数据
tmp=zeros（1,M）;
%初始化中间变量，用于保存每段观测数据的功率
fori=1:
L%观每段有M=N/L个数据
forj=1:
M
y（i,j）=x（j+i*M-M）;
%取出L段数据
end
xk=fft（y（i,:
）,M）;
%fft变换
Ix=（（abs（xk））.^2）./M;
%计算每段的功率
Px=tmp+Ix;
%累加每段的功率
end
Px_L=Px./L;
%平均功率谱
二、实验题目二
（1）、信号产生
首先根据题目给定的信噪比计算正弦信号的幅度，产生所需的正弦信号后与噪声信号叠加形成观测信号。
若正弦信号幅度为A，则其功率为：
信噪比：
，其中P为噪声功率；
由以上公式便可计算出两个正弦信号的幅度。
（2）利用MATLAB功率谱分析
自相关法
首先由观测数据估计出其自相关函数，再解该方程得到模型参数，进而求得信号的功率谱，一般利用Levenson-Durbin递推法求解。
matlab中的自相关法谱估计如下：
pyulear（xn,p,nfft,fs）;
其中，xn是输入观测信号，p是阶数，nfft是观测数据两倍长度的频率采样点，fs是采样频率。
Burg法
Burg法直接利用观测数据计算AR模型参数，适用于观测数据长度较短的时候。
matlab中Burg法谱估计如下：
pburg（xn,p,nfft,fs）;
其中，px1是返回的功率谱，f1是返回的频率分布；
xn是输入观测信号，p是阶数，nfft是观测数据两倍长度的频率采样点，fs是采样频率。
协方差法
协方差法利用到的数据完全一致，相对自相关法没有数据为0的假设。
matlab中协方差法谱估计如下：
pcov（xn,p,nfft,fs）;
修正协方差法
修正协方差法使用到了前向和后向预测误差对模型参数进行估计。
matlab中的修正协方差法谱估计如下：
pmcov（xn,p,nfft,fs）;
（1）数据长度的影响：
将四种方法产生的功率谱显示在了一幅图中。
L=1024时，如下图9所示。
L=256时，结果如下图10所示。
图9L=1024
图10L=256
分析：
对比两幅图可以发现：
随着数据长度的减小功率谱的分辨率在逐渐降低。
即数据长度越大，功率谱分辨率越高，数据长度越小，功率谱分辨率越低。
（2）模型阶次的影响：
当阶次发生变化时，如p=100时结果如图11所示，p=30时，结果如图12所示。
图11阶次p=100
图12阶次p=30
对比阶次p=100与阶次p=30时的结果可以看出，当阶次太小时，分辨率比较低，无法估计出小信噪比的信号，即小信噪比信号无法在图中观察出来。
（3）信噪比的影响：
将110Hz正弦信号信噪比降低到10，结果如下图13所示：
图13SNR1=10SNR2=10
对比发现：
变化最明显的是自相关法和Burg法，当两个正弦信号信噪比相差不多时，自相关法的估计效果还是较好的；
而Burg法效果很差，每次实验的结果不同，经常分辨不出一个正弦频率，出现了谱峰偏移和谱线分裂的现象。
而对于自相关法，并不适合于分辨信噪比相差较多的两个正弦信号混合的情况，结果会谱峰偏移，出现观测不到一个正弦信号的频率。
如修改信噪比SNR1=30SNR2=10的情况，结果如下图14所示，仍然分辨不出信噪比较低的信号，即110Hz，SNR2=10的正弦信号。
图13SNR1=30SNR2=10
1、通过实验，加深了对AR模型谱估计的流程，掌握了结合matlab利用自相关法、Burg法、协方差法、修正协方差法估计功率谱，体会了各种方法的性能，优缺点等；
2、能更加熟练地通过matlab帮助实现对函数的查找，使用；
3、由于考试最近考试较多，时间有限，没有自己编程实现各种现代功率谱估计的方法，在时间空闲时，也会结合课堂讲授的流程图对各种方法尝试自行编程实现。
附matlab源代码
实验题目一
clc;
%清除命令窗口
clear;
%清空变量
N=1024;
%数据长度
df=1;
%设定频率分辨率
fs=（N-1）*df;
%计算采样率
t=0:
1/fs:
（N-1）/fs;
%截取信号的时间段
f=0:
df:
fs;
%功率谱估计的频率分辨率和范围
%平稳随机信号的获取
%准确获得观测信号是很关键的
%处理方法一：
迭代的方式
x（i）=0.8*x（i-1）+w（i）;
%处理方法二：
卷积的方式
%n=1:
%h（n）=（0.8）.^（n-1）;
%要注意n-1不是n，因为冲击响应是从0开始的
%y=conv（w（n）,h（n））%共2*N-1项，取前N项即可
%令N=5时，实验结果，取卷积的前N项，可以看出两种方式结果是相同的
%x=
%0.62321.29761.97900.59110.6849
%y=
%0.62321.29761.97900.59110.68490.34370.0131
%-0.29780.0868
%当数据长度为256时，取前256个数据为x2
N2=256;
x2=x（1:
256）;
%理想功率谱计算
H=fft（h,N）;
%对冲击响应做傅里叶变换
Px_ideal=（abs（H））.^2;
Px_ideal_db=10*log10（Px_ideal）;
%以分贝为单位表征
%周期图法估计功率谱
%数据长度为1024时
xk=fft（x,N）;
%傅里叶变换
Px_zqt=（（abs（xk）.^2）./N）;
Px_zqt_db=10*log10（Px_zqt）;
%数据长度为256时
xk2=fft（x2,N）;
Px_zqt2=（（abs（xk2）.^2）./N2）;
Px_zqt_db2=10*log10（Px_zqt2）;
%Px_mean（L,x）:
平均周期图法,分段数目L，数据为x
Px_1024_2=Px_mean（2,x）;
%参数L=2x（1024个数据）时
Px_1024_8=Px_mean（8,x）;
%参数L=8x（1024个数据）时
Px_256_2=Px_mean（2,x2）;
%参数L=2x（256个数据）时
Px_256_8=Px_mean（8,x2）;
%参数L=8x（256个数据）时
Px_1024_2_db=10*log10（Px_mean（2,x））;
Px_1024_8_db=10*log10（Px_mean（8,x））;
Px_256_2_db=10*log10（Px_mean（2,x2））;
Px_256_8_db=10*log10（Px_mean（8,x2））;
%显示噪声w（n）与平稳信号x（n）
figure;
plot（n,w（n））;
xlabel（'
n'
）;
ylabel（'
w（n）'
legend（'
随机噪声w（n）'
）
plot（n,x（n））;
x（n）'
平稳随机信号x（n）'
%显示理想功率谱
plot（f,Px_ideal）;
xlabel（'
频率f（Hz）'
ylabel（'
功率谱Px_ideal'
理想功率谱Px_ideal'
plot（f,Px_ideal_db）;
功率谱Px_ideal_db（dB）'
理想功率谱Px_ideal_db'
%周期图法估计功率谱结果
subplot（2,2,1）;
plot（f,Px_zqt）;
功率谱Px'
1024'
subplot（2,2,2）;
plot（f,Px_zqt_db）;
功率谱Px（dB）'
subplot（2,2,3）;
plot（f,Px_zqt2）;
256'
subplot（2,2,4）;
plot（f,Px_zqt_db2）;
%平均周期图法结果
plot（1:
length（Px_1024_2）,Px_1024_2）;
L=2N=1024'
length（Px_1024_8）,Px_1024_8）;
L=8N=1024'
length（Px_256_2）,Px_256_2）;
L=2N=256'
length（Px_256_8）,Px_256_8）;
L=8N=256'
子程序Px_mean.m
%计算数据长度
实验题目二
fs=1000;
%采样频率
n=0:
2;
%取样数据点
snr1=30;
%正弦波d的信噪比
snr2=10;
%正弦波d2的信噪比dB
A1=sqrt（2*10^（snr1/10））;
%正弦波x1的幅度
A2=sqrt（2*10^（snr2/10））;
%正弦波x2的幅度
An=sqrt
（1）;
%噪声w（n）的幅度
x1=A1*sin（2*pi*100*n）;
%正弦信号x1
x2=A2*sin（2*pi*110*n）;
%正弦信号x2
w=An*randn（size（n））;
%噪声信号w（n）
xn=x1+x2+w;
%观测信号
L=1024;
%数据长度
xn=xn（1:
L）;
%取L长度的观测数据
n=n（1:
size（xn,2））;
%取得数据点
nfft=2*L;
%频域数据长度
p=80;
%模型阶次
%自相关法功率谱估计
[px1,f1]=pyulear（xn,p,nfft,fs）;
plot（f1,px1）;
title（'
自相关法功率谱估计'
axis（[8015001.2*max（px1）]）;
%协方差法功率谱估计
[px2,f2]=pcov（xn,p,nfft,fs）;
plot（f2,px2）;
协方差法功率谱估计'
axis（[8015001.2*max（px2）]）;
%修正协方差法功率谱估计
[px3,f3]=pmcov（xn,p,nfft,fs）;
plot（f3,px3）;
修正协方差法功率谱估计'
axis（[015001.2*max（px3）]）;
%burg法功率谱估计
[px4,f4]=pburg（xn,p,nfft,fs）;
plot（f4,px4）;
burg法功率谱估计'
axis（[8015001.2*max（px4）]）;
升级会员 