数字信号处理之信号处理仿真 第五章 功率谱估计(2)

之前讲解了周期图法，也给出了相关仿真程序！由浅入深的学！但是大家可能已经发现了一个现象：周期图法估计出的功率谱不够精细，分辨率比较低。因此需要对周期图法进行修正，可以将信号序列x(n)分为n个不相重叠的小段，分别用周期图法进行谱估计，然后将这n段数据估计的结果的平均值作为整段数据功率谱估计的结果。还可以将信号序列x(n)重叠分段，分别计算功率谱，再计算平均值作为整段数据的功率谱估计。这两种方法称为分段平均周期图法，一般后者比前者效果好。加窗平均周期图法是对分段平均周期图法的改进，即在数据分段后，对每段数据加一个非矩形窗进行预处理，然后在按分段平均周期图法估计功率谱。相对于分段平均周期图法，加窗平均周期图法可以减小频率泄漏，增加峰值频率的宽度。后面看看分段平均周期图法的仿真程序，看看谱估计的效果！预先介绍一下pwelch函数。初学者在学习谱估计知识的时候可以借助MATLAB自带函数的力量实现快速成长！当年我就从中收益颇多！了解函数一般是从MATLAB的help入手，当然很多人也愿意从百度入手，效果一样！

pwelch函数利用Welch平均功率图法返回信号x的功率谱密度（PSD）。

当x是向量时，它被当做一个单通道信号。当x是矩阵，x的每一列被当做一个通道的信号，其psd结果相对应与psd的每一列。

如果x是实值信号，则pxx是单边谱估计，如果x是复数信号，则pxx是双边谱估计。

在默认设置中，x被分成8段，重叠率为50%的片段。每个片段上添加一个hamming窗。pwelch利用这种修正后的周期图法来对psd进行估计。

如果不能将x刚好分为满足50%重叠的8个片段，为了实现功能，函数会对x的长队进行相应地自动裁剪。

pwelch函数中参数的具体意义：

x：进行功率谱估计的有限长输入序列；WINDOW：指定窗函数，默认值为hamming窗；NFFT：DFT的点数；Fs ：绘制功率谱曲线的抽样频率；Pxx：功率谱估计值；F：Pxx值所对应的频率点；NOVERLAP指定分段重叠的样本数，如果NOVERLAP=L/2，则可得到重叠50%的Welch法平均周期图。如果使用boxcar窗，且NOVERLAP=0，则可得到Bartlett法的平均周期图。

[pxx,w] = pwelch(___)返回的是标准化后的频率向量w，如果pxx是单边谱估计，那么w的范围就是0到pi，如果pxx是双边谱估计，那么w的范围就是0到2pi。

[pxx,f] = pwelch(___,fs) 返回一个频率向量f。fs时每个单位时间的样本s，那么f的单位是Hz，那么f的范围是0到fs/2，对于复值信号来说，f的范围是0到fs。且fs必须是pwelch的第五个输入变量，如果输入了fs，而其他参数使用默认值的话，其他参数可以设置为空。

先看之前文章提出的问题！如果信号的频率很接近，谱估计会带来什么结果？这就涉及了分辨率的概念！

似乎间接法分辨的效果也不是很好了！

那怎么办呢？增加数据的采样时间是最好的也是最直接的办法！和FFT需要增加信号时长的道理是一致的。好了，该看代码了！本文的程序非常的长，希望大家有耐心的看完！因为需要进行多种谱估计方式的结果比对，自然程序内容就多了起来！

% psd_sim2.m

% 采用经典谱估计中的周期图法分析信号的功率谱！

% 采用直接法与间接法两种方式产生信号的功率谱。

% 测试频率相近信号下的谱估计结果！

% 软件版本：2021a

clear all;close all;

%%%%%%%%%  生成信号

Fs = 1000;  % 采样频率

%%% 产生数字信号

time = 0:1/Fs:1-1/Fs; 

% 时间序列

f1 = 50;   % 信号频率

f2 = 55;   % 信号频率

% 注明：采样频率的设置必须符合奈奎斯特准则！

signal = cos(2*pi*f1*time) + 3*cos(2*pi*f2*time) + randn(1,length(time));

subplot(5,1,1);

plot(signal);

title('加噪信号');

grid on

%%%%%%%%%  周期图法