现代数字信号处理实验

现代信号处理作业 实验题目实验题目 设信号 8 0cos 25 0 47 0 cos 35 0 cos 321 nvnnnnx 其中是内的独立随机变量 v n 是单位高斯白噪声 321 1 利用周期图法对序列进行功率谱估计 数据窗采用汉明窗 2 利用 BT 法对序列进行功率谱估计 自相关函数的最大相关长 度为 M 64 128 256 512 采用 BARTLETT 窗 3 利用 Welch 法对序列进行功率谱估计 50 重叠 采用汉明 窗 L 256 128 64 4 利用 Burg 法对序列进行 AR 模型功率谱估计 阶数分别为 10 13 要求每个实验都取 1024 个点 fft 作为谱估计 取 50 个样本序 列的算术平均 画出平均的功率谱图 实验原理 1 周期图法 又称间接法 它把随机信号的 N 个观察值 xN n 直接进行傅 里叶变换 得到 XN ejw 然后取其幅值的平方 再除以 N 作为对 x n 真实功率谱的估计 其中 2 1 jw eX N wP N per 1 0 N n jwn N jw N enxeX 2 BT 法 对于 N 个观察值 x 0 x 1 x N 1 令 xN n a n x n 计算 rx m 为 计算其傅里叶变换 mN n NNx Nmmnxnx N mr 1 0 1 1 作为观察值的功率谱的估计 M Mm jwm xBTNMemrmvwP1 其中 v m 是平滑窗 3 Welch 法 假定观察数据是 x n n 0 1 2 N 1 现将其分段 每段长度为 M 段与 段之间的重叠为 M K 第 i 个数据段经加窗后可表示为 1 1 0 MiiKnxnanxiM 其中 K 为一整数 L 为分段数 该数据段的周期图为 其中 由此得到平均周期图为 2 1 wX MU wP i M i per 1 0 M n jwni M i M enxwX 其中归一化 U 取 1 0 1 L i i perwP L wP 1 0 2 1 M n na M U 4 Burg 法 在约束条件下 使得极小化 其中 约束条件是它所得到 2 1 bf 的各阶模型解要求满足 Levison 递归关系 仿真结果仿真结果 1 周期图法 2 1 BT 法 平滑窗采用 BARTLETT 窗 长度为 64 2 2 BT 法 平滑窗采用 BARTLETT 窗 长度为 128 2 3 BT 法 平滑窗采用 BARTLETT 窗 长度为 256 2 4 BT 法 平滑窗采用 BARTLETT 窗 长度为 512 3 Welch 法 3 1 L 256 3 2 L 128 3 1 L 64 4 1 Burg 法 阶数为 10 4 1 Burg 法 阶数为 13 仿真分析仿真分析 1 周期图法得到的功率谱 特点是离散性大 曲线粗糙 方 差较大 但是分辨率较高 2 从图中可以看出间接法估计的偏差大于周期图法 这是因 为 BT 法在上施加了一个较短的平滑窗 v m 但是 BT 法得到 mrx 的功率谱的方差 从图中可以看出 要小于周期图法得到的方差 所以其分辨率要比周期图法差 从理论上 BT 法和周期图法的方差之比为 一般 M Mn nv N 1 2 来说 v m 是以 m 0 对称递减的 又 M N 所以 1 验证了实验 所得 BT 法中对于延迟窗取不同的长度 从图中可以看出 长度越长 方差越大 分辨率越高 因为 BT 法的方差和成正比关系 M Mn nv 2 当长度越长时 越大 所以得到的功率谱的方差越大 M Mn nv 2 BT 法中要求平滑窗的长度为 2M 1 即为奇数 关于 m 0 处 对称 本实验中 给出的平滑窗为偶数 与自相关函数关于 m 0 对 称 且 r 0 最大这一特性是否十分吻合 值得商榷 3 Welch 法是将 N 点观察值分为 L 个数据段 用以改善功率 谱图的方差特性 从图中看出 对比周期图法 Welch 法作出的功 率谱图的方差特性的确得到大大的改善 但是在给每段序列用适当 的窗口函数加权后 在得到平滑的估计结果同时 使得功率谱额主 瓣变宽 因此分辨率有所下降 从图中可以看出 Welch 法中 随着分的段数增加 得到的功 率谱的方差变得更好 这是因为 Welch 法中方差与分的段数大约成 反比关系 而分辨率则随着段数增加而下降 4 burg 法从图中可以看出 随着阶数的增高 分辨率提升 本实验中 由于阶数比较小 得到的功率谱图并不理想 现代谱估 计的一些隐含着数据和自相关函数的外推 使其可能的长度超过给 定的长度 不象经典谱估计那样受窗函数的影响 因而现代谱的分 别率比较高 而且现代谱线要平滑得多 附录 附录 1 周期图法 N 512 Nfft 1024 n 0 N 1 N1 50 e1 unifrnd pi pi e2 unifrnd pi pi e3 unifrnd pi pi w hamming N yn zeros Nfft N1 zn zeros Nfft 1 for i 1 N1 xn cos 0 35 pi n e1 cos 0 4 pi n e2 0 25 cos 0 8 pi n e3 randn 1 N px abs fft w xn Nfft 2 1024 yn i px end for i 1 N1 zn zn yn i end znf 10 log10 zn N1 f 0 length px 1 length px plot f znf xlabel frequency hz ylabel power spectrum Db title periodogram N 512 2 bt 法 N 512 M 64 平滑窗的长度 可更改 Nfft 1024 n 0 N 1 N1 50 e1 unifrnd pi pi e2 unifrnd pi pi e3 unifrnd pi pi w1 bartlett M w2 zeros 1 512 M 2 length w2 w3 zeros 1 Nfft 1 M w w2 w1 w2 生成 BARTLETT 窗 length w yn zeros Nfft N1 zn zeros Nfft 1 for i 1 N1 xn cos 0 35 pi n e1 cos 0 4 pi n e2 0 25 cos 0 8 pi n e3 randn 1 N rx xcorr xn biased 求序列自相关 rx rx 0 px abs fft w rx Nfft 自相关函数与平滑窗做 fft yn i px end for i 1 N1 zn zn yn i end znf 10 log10 zn N1 f 0 length px 1 length px plot f znf xlabel frequency hz ylabel power spectrum Db title BT N 64 3 welch 法 N 512 Nfft 1024 n 0 N 1 N1 50 e1 unifrnd pi pi e2 unifrnd pi pi e3 unifrnd pi pi L 3 生成 15 个数据段 M 2 N L 1 每段数据的长度 w hamming M U w w M 归一化因子 yn zeros Nfft N1 zn zeros Nfft 1 px zeros Nfft 1 for i 1 N1 xn cos 0 35 pi n e1 cos 0 4 pi n e2 0 25 cos 0 8 pi n e3 randn 1 N for i 1 L px1 abs fft w xn i 1 2 N 2 L 2 1 i 1 2 N 2 L 2 M Nfft 2 M U 每一段进行周期图 px px px1 end px2 px L 得到平均周期图 yn i px2 end for i 1 N1 zn zn yn i end f 0 length px 1 length px znf 10 log10 zn N1 plot f znf xlabel frequency hz ylabel power spectrum Db title Welch L 256 4 burg 法 Nfft 1024 n 0 N 1 N1 50 e1 unifrnd pi pi e2 unifrnd pi pi e3 unifrnd pi pi yn zeros Nfft 2 1 N1 length yn zn zeros Nfft 2 1 1 for i 1 N1 xn cos 0 35 pi n e1 cos 0 4 pi n e2 0 25 cos 0 8 pi n e3 randn 1 N order 10 range half px1 w pburg x