这两天在做一个《数字信号处理》的作业，主要是功率谱估计方面的知识运用，这里又不得不再次检讨一下自己的磨蹭了，一个期末大作业居然做了整整两天，简直慢到怀疑人生。
话不多说，上题目。
一、试验数据的产生
分别产生两个零均值的高斯白噪声数据u1(n)和u2(n)，其功率都为 σ2 = 0.12，让 u1(n)和u2(n)分别通过一个FIR系统，得到输出为v1(n)和v2(n)。该FIR系统由 5个FIR子系统级联而成：
B1 = [1;1.98;0.98]，
B2 = [1; −1.98;0.98]，
B3 = [1; −1.8418;0.98]，
B4 = [1; −1.5;0.98]，
B5 = [1; −1.2727;0.98]，
令v(n) = v1(n) + jv2(n)。
在v(n)上再加上三个复正弦信号，它们的幅度分别为a1=12;a2 =12;a3 =6，归一化频率分别为f1 = 0.22; f2 = 0.21; f3 = 0.12，这样可得到已知功率谱的试验信号x(n) 。
1. 写出得到 x(n) 的计算过程。令所得到的数据长度 N = 256，描绘该波形。
2. 描绘出该试验信号的真实功率谱Px(ejw)。Px(ejw)在一个周期内取4096个点（注：以下凡是功率谱曲线，都这样取）。
3. 求x(n)中三个正弦信号在上述三个频率处的信噪比（注：为简化考虑，此处信噪比仅考虑正弦对白噪声，而不考虑白噪声通过系统后的影响）。
二、 对x(n)做功率谱估计
1 用周期图估计 x(n) 的功率谱，令N = 256，画出Pper(ejw)曲线。
2. 用 Welch 方法对上述周期图法作改进，例如，可选每段64点，重叠32点，用 Hamming 窗，画出所求功率谱曲线。
3. 用自相关法估计x(n)的功率谱，令 M = 63，画出PBT(ejw) 曲线
4. 用 AR 模型估计 x(n) 的功率谱，阶次p可在6到15之间，由自己给定，给出估计的功率谱曲线 PAR( ejw) 。分别用自相关法和Berg算法。
5. 试利用 FPE 或 AIC 准则来确定本信号所需要的 AR 模型的阶次。
题解：
1、白噪声信号的产生
利用rand产生一组白噪声数据。题目中要求的白噪声数据长度是N=296。为了保证产生的数据能较好的近似白噪声，数据的长度要尽量取得长一些，此处取原始数据长度NN=2000。将得到的白噪声数据先减去均值，再乘以 σ ，即可得到满足要求的白噪声数据 。同理，再次利用上述方法可以生成另一组白噪声 ，因为是计算机随机产生的两组数据，所以两组序列不相关如果程序中取 σ2 =0.12。
2、系统级联
利用MATLAB的conv函数，可以求得级联系统为11阶FIR系统，系统响应为
归一化系统响应为：
3、试验信号的叠加
令u1(n)和u2(n) 分别通过该级联系统，分别输出v1(n) ，v2(n) 。注意，在计算u(n) 与h(n)的卷积时可以调用零相滤波器函数filtfilt避免时延对系统的影响。
再令
v(n) = v1(n) + jv2(n)
则 是一个复值的噪声序列，其实部和虚部也基本不相关，那么 的真实功率谱可以求出，即
其中b(k)是模型的系数。
在v(n)的基础上加上三个复正弦，其归一化频率已知。f1 = 0.22; f2 = 0.21; f3 = 0.12。即
然后可以取中间的数据段截取N=256个点这样就可以求出试验信号xn了。
4、求出xn的真实功率谱。复正弦信号的强度是piA2 。因此，xn的真实功率谱是
下图所示即为xn的真实功率谱，其归一化频率范围是从-0.5~0.5，
5、求x(n)中三个正弦信号在上述三个频率处的信噪比。根据信噪比的定义
其计算方法是10log(Ps/Pn)，其中Ps和Pn分别代表信号和噪声的有效功率。
那么：
将数据代入公式得
第一个正弦信号的信噪比SNR1 = 27.7815
第二个正弦信号的信噪比SNR2 = 27.7815
第三个正弦信号的信噪比SNR3 = 21.7609
6、用周期图法对试验信号进行谱估计。其原理是把随机信号xn的N点观察数据视为能量有限信号，直接进行傅立叶变换，然后再取其幅值的平方，并除以N，具体公式如下
在MATLAB中可以直接调用periodogram函数来计算功率谱密度。
如图所示，由于主瓣宽度过小，f1和f2不能完全分开，只是在波形的顶部能看出是两个频率的分量。
7、用自相关法对试验信号进行谱估计。其原理是先估计出 xN(n)自相关函数r(m)，然后对r(m)求傅立叶变换得到xN(n)的功率谱，以此作为P(w)的估计。
下图是自相关法估计的效果图，显然谱变得光滑，但是较大的边瓣使分辨率降低了很多。
8、用Welch法估计试验信号的功率谱，Welch方法是一种修正周期图功率谱密度估计方法，它通过选取的窗口对数据进行加窗处理，先 xN(n)分段时允许数据重叠求功率谱之后再进行平均。
MATLAB中可以直接调用pwelch函数估计信号的功率谱，其基本的调用格式
[Pxx,f]=pwelch(x,window,noverlap,nfft,fs)
其中窗口的长度N表示每次处理的分段数据长度，Noverlap是指相邻两段数据之间的重叠部分长度。
根据题目要求，选每段64点，重叠32点，用 Hamming 窗，画出所求功率谱曲线如下图所示。谱变得更加平滑，但是分辨率相比周期图法有所降低。
9、用AR模型估计x(n)的功率谱。其最常用的有两种方法，自相关法和Burg法。自相关法的基本思想是利用Levinson递推公式求解Yule-walker公式求得的AR模型的参数等效于前向预测器的参数。因此AR模型的自相关法等效于对前向预测序列efp(n)前后加窗；Burg算法的基本思想是直接从观测的数据利用线性预测器的前向和后向预测的总均方误差之和为最小的准则来估计反射系数，进而通过Levinson递推公式求出 AR模型优化参数。
根据本题的数据可以得出下面的AR模型估计谱，其中自相关法估计的阶数取的是15阶，Burg算法取得是6阶，这里可以看出无论是在谱的分辨率方面还是模型的阶次方面，Burg算法都比自相关算法有明显的优点。经过试验，自相关算法的阶次大概取30阶时才能把 和 两个谱峰分开，但是Burg算法在阶次等于5的时候就能完全把两个谱峰完全分开。
10、利用FPE准则来确定本信号所需要的AR模型的阶次。因为从低阶到高阶，阶次选得高，模型的最小预测误差功率是递减的，固然会提高谱估计的分辨率，但同时会产生虚假的谱峰或谱的细节。 所以必须选择合适的AR模型的阶次。这里可以根据最终预测误差准则(FPE)来确定AR模型的阶次。当阶次k由1增加时，FPE(k)将在某一个k处取一个极小值，可以将此时的k定为最合适的阶次p。具体的计算公式如下：
在MATLAB中可以调用aryule函数和arburg函数来确定k阶预测器的最小预测误差 ，实际计算
Optimal_FPE_yule =16；Optimal_FPE_burg=49；
但是可以画出FPE(k)的取值趋势变化曲线，可以看到，实际上使用Burg算法时，在阶次大于5后 的变化不大，所以实际计算中可以直接取6阶计算。这也与上一步的试验分析一致。
