功率谱估计及比较

-60 -80 -100 -120 -140 -160 -180 -200 -0.5
-0.4
-0.3
-0.2
-0.1 0 0.1 归一化频率序号k
0.2
0.3
0.4
0.5
图 1 周期图法求功率谱（左：编程 右：函数）
可以看出两种方法的计算结果相同，也验证了编程的有效性。同时，对应于 周期图法所求得的功率谱振荡剧烈，信号方差较大，不利于对功率信号的分析。 (2) 自相关法 先根据实验所给数据求解出自相关函数， 然后对自相关函数进行傅里叶变换， 从而得到功率谱估计。自相关法是由维纳-辛钦公式出发的，本质上是对周期图 法的插值，因此而这本质上来说是一致的。从图 2 可以看出，自相关法得到的结 果与周期图法相似，同样是方差较大，信号振荡大。
(3) 不同的窗函数对 Welch 谱估计结果有何影响？
归 一 化 功 率 谱 信 号 S /dB
不 同 窗 函 数 下 的 Welch谱 估 计 函 数 ， 叠 合 长 度 为 50% ， 每 段 长 度 为 100 0 矩形窗 -10 三角窗 布莱克曼窗 汉明窗 -20 汉宁窗 -30 -40 -50 -60 -70 -80 -90 -0.5
功率谱估计分析及比较
1 实验目的
(1) 掌握Welch算法的概念、应用及特点； (2) 了解谱估计在信号分析中的作用； (3) 能够利用 Welch 法对信号作谱估计，对信号的特点加以分析。
2 实验内容
(1) 读入实验数据。 (2) 编写一利用Welch法作估计的算法程序。 (3) 将计算结果表示成图形的形式，给出信号谱的分布情况图。
4 实验结果分析
4.1 实验仿真结果 在实验仿真求取结果中， 可以采取两种方式， 一是采用第 3 部分的算法步骤， 进行编程求解；二是在 Matlab 软件中调用 periodogram.m、pwelch.m 等已有的函 数进行计算求解。 (1) 周期图法 函数介绍：[Sxx,f]=periodogram(x,window,nfft,fs,'range') x 为进行功率谱估计的输入有限长序列，window 用于指定采用的窗函数， nfft 设定 fft 算法的长度，fs 为采样频率，Sxx 为输出的功率谱密度值，f 为得到 的频率点。 分别采用编程和函数两种方法求解所有数据的功率谱， 并且进行归一化后的 结果如下图所示：
I i ( )

1 M 1 | x(n)e jm |2 M n 0
i 1, 2,..., k
k 段平均后的谱估计为：
S xx ( )

1 k I i () k i 1
(5) 韦尔奇（Welch）谱估计法 韦尔奇谱估计法将加窗平滑法和平均周期法二者相结合，其原理在于： a) 把 N 个数据分成 k 段，每段可以互相独立(如平均周期图法)，也可以互 相交叠，例如交叠一半，即 k N / L 或 k ( N L / 2) / ( L / 2) ； b) 再把每段数据乘上窗函数 w(n) (如加窗平滑法)后作 DFT。 其操作步骤如下： a) 每段数据长度 L 补 0 至 M（2 的整数幂）个点； b) 第 i 段 xi (n) 的长度 L M
对于直接法的功率谱估计，当数据长度 N 太大时，谱曲线起伏加剧，若 N 太小，谱的分辨率又不好，Welch 谱估计法从两个方面对直接法进行改进： 一是分段，且各段之间可以存在区间重叠，这样会使方差减小； 二是加窗， 选择适当的窗函数加在 x(n) 的分段数据上， 这样可使谱估计变得 更加平滑，但分辨率减小。 从理论分析上来说， Welch 谱估计法是渐进无偏估计， 并且是渐进一致估计。
3 算法讨论与分析
随机信号的功率密度函数的定义如下：
S xx ( j) lim E[| X T ( j) |2 ]
T
其物理意义表示为信号在单位频带内的平均功率。 经典的谱估计主要包括以 下几种方式： (1) 周期图法 周期图法是直接建立在功率谱的定义式上的，也称之为直接法。原理计算如 下：
归 一 化 功 率 谱 信 号 S /dB
-80 -100 -120 -140 -160 -180 -200 -0.5
-0.4
-0.3
-0.2
-0.1 0 0.1 归一化频率序号k
0.2
0.3
0.4
0.5
5 原程序清单
clc; clear; Qjt=load('Qjt.dat'); %读取数据 n=length(Qjt); nfft=2^(nextpow2(n)-1); fs=50; data(:,1)=Qjt(1:nfft,1); %---------------周期图法（编程求解）------------y=fft(data); Sxx=abs(y).^2/n; Sxx=Sxx/max(Sxx); S_per2=fftshift(20*log10(Sxx)); f=-0.5:1/nfft:0.5-1/nfft; figure(12) plot(f,S_per2,'b'); title('直接法（编程）求功率谱密度函数') xlabel('归一化频率序号 k'); ylabel('归一化功率谱信号 S-per /dB') grid; %----------周期图法（函数：periodogram）--------[Sxx]=periodogram(data,[],nfft,1,'twosided'); %给定参数下对数据进行双边周期图法功率谱估计 Sxx=Sxx/max(Sxx); S_per1=fftshift(20*log10(Sxx)); f=-0.5:1/nfft:0.5-1/nfft; figure(11) plot(f,S_per1,'b') title('直接法（函数）求功率谱密度函数') xlabel('归一化频率序号 k'); ylabel('归一化功率谱信号 S-per /dB') grid; %------------------自相关法--------------------rxx=xcorr(data,'biased'); Sxx=fft(rxx); Sxx=Sxx/max(Sxx); S_bt=fftshift(20*log10(Sxx)); f=-0.5:1/2/nfft:0.5-1/nfft; figure(2) plot(f,S_bt,'b') title('自相关法求功率谱密度函数'); %求解自相关函数 % 做 FFT，求解功率谱 % 将功率谱进行归一化 %数据化为分贝值，并进行频移 %归一化频率 % 将功率谱进行归一化 %数据化为分贝值，并进行频移 %归一化频率 % 作图 % 对原始数据做 FFT % 求模平方并除以 N，即为功率谱 %将功率谱进行归一化 %数据化为分贝值，并进行频移 %归一化频率 % 作图 % 给定数据的长度 %截取最大的 2 的整次幂数据长度 % 信号数据采样频率 %为加快计算，选用数据长度为 2 的整次幂
该方法的计算步骤： a) 取 N 点数据的 DTFT（DFT） b) 求模之平方并除以 N (2) 自相关法 自相关法的原理是由维纳-辛钦公式，经自相关函数间接获得的。原理计算 如下：
1 N |m|1 R ( m ) xx x ( n ) x ( n m) N n 0 2 N 1 j km S BT (k ) N R ( m ) e xx m ( N 1)
1 a) 作 R xx (m) N
N |m|1


n 0
x(n) x(n m)
m 0, 1, 2,..., M
b) 选取 v(m) c) 作 S BT ( )

m M
v(m) R
M

xx
(m)e jm
(4) 平均周期图法（Bartlett 法） 平均周期图法的原理是 k 个独立同分布的随机变量的均值之方差， 等于单个 变量方差的 1/ k 。 具体的方法步骤是长数据 N 分成 k 段，每段 M N / k ，针对每段分别用周 期图法求谱：
归 一 化 功 率 谱 信 号 S-wel /dB
-0.4
-0.3
-0.2
-0.1 0 0.1 归一化频率序号k
0.2
0.3
0.4
0.5
4.2 结果分析 (1) 不同分段点数对 Welch 谱估计结果有何影响？ 重叠长度为 50%，窗函数为矩形窗，分别取分点点数为 100、500、1000。
不 同 分 段 点 数 的 Welch谱 估 计 函 数 ， 叠 合 长 度 为 50% 0 -10 -20 L=100 L=500 L=1000
自相关法求功率谱密度函数 0 -20 -40
归 一 化 功 率 谱 信 号 S-bt /dB
-60 -80 -100 -120 -140 -160 -180 -200 -0.5
-0.4
-0.3
-0.2
-0.1 0 0.1 归一化频率序号k
0.2
0.3
0.4
0.5
图 2 自相关法求功率谱函数
(3) 加窗平滑法（BT） 根据加窗平滑法的求解步骤进行编程，取窗函数为矩形窗，其长度为 M。对 应于不同的 M 值可以得到不同的结果，如图三所示。
-0.4
-0.3
-0.2
-0.1 0 0.1 归一化频率序号k
0.2
0.3
0.4
0.5
(4) 周期图法与 Welch 法的谱估计结果有何不同？
周 期 图 法 和 Welch法 谱 估 计 对 比 （ 每 段 长 度 100， 叠 合 50% ， 矩 形 窗 函 数 ） 0 S-per -20 S-welch -40
-10
归 一 化 功 率 谱 信 号 S /dB
-20
-30
-40
-50
-60
-70 -0.5
-0.4
-0.3
-0.2
-0.1 0 0.1 归一化频率序号k
0.2
0.3
0.4
0.5
可以分析得到， 如果分段的数据互补交叠， 则估计方差和平均周期图法一样。 如果数据相互交叠，则分段数 K 比平均周期图法增大，则估计方差比平均 周期图法小。 实验仿真结果也可以看出，交叠长度仅为 1 是其振荡叫交叠长度为 69 的略 微剧烈，故其估计方差较大。 值得注意的是，数据的相互交叠，也降低了数据段之间的独立性，实际方差 的减少也不如理想中那么好。
平均周期图法求功率谱密度函数 0 -10 -20
归 一 化 功 率 谱 信 号 S-bar /dB
-30 -40 -50 -60 -70 -80 -90 -0.5
-0.4
-0.3
-0.2
-0.1 0 0.1 归一化频率序号k
0.2
0.3
0.4
0.5
图 4 平均周期图法求功率谱（K=128）
(5) Welch 法 函数介绍： [Px,F]=pwelch(x,window,Noverlap,Nfft,fs) 分别采用编程和函数两种方法求解所有数据的功率谱， 并且进行归一化后的 结果如下图所示：
加窗平滑法求功率谱 0 -20 -40 M=50 M=100 M=1000
归 一 化 功 率 谱 信 号 S /dB
-60 -80 -100 -120 -140 -160 -0.5
-0.4
-0.3
-0.2
-0.1 0 0.1 归一化频率序号k
0.2
0.3
0.4
0.5
图 3 加窗平滑法求功率谱函数
(4) 平均周期图法（Bartlett） 将数据平均分为 K 段
该方法的计算步骤： a)
x(n) ，N 点 2 N 1点，得 R xx (m)

b) 按 2N-1 点对 R xx (m) 作 DFT， R xx (m) S BT (k ) (3) 加窗平滑法（BT 法） 加窗平滑法的原理是先做自相关估计， 在选择合适的窗函数相乘， 也即截断， 然后作 DFT。其原理步骤如下：
归 一 化 功 率 谱 信 号 /dB
-60 -80 -100 -120 -140 -160 -180 -200 -0.5
-0.4
-0.3
-0.2
-0.1 0 0.1 归一化频率序号k
0.2
0.3
0.4
0.5
所有谱估计法所求解的功率谱 （ 加 窗 平 滑 窗 函 数 长 度 为 50， 平 均 周 期 图 法 均 分 为 128段 Welch法 每 段 数 据 长 度 100， 叠 合 50% ， 矩 形 窗 ） 0 -20 -40 -60 周期图法 自相关法 加窗平滑法 平均周期图法 Welch法
直接法（编程）求功率谱密度函数 0 -20 -40 0 -20 -40 直接法（函数）求功率谱密度函数
归 一 化 功 率 谱 信 号 S-per /dB
-60 -80 -100 -120 -140 -160 -180 -200 -0.5
归 一 化 功 率 谱 信 号 S-per /dB
-0.4 -0.3 -0.2 -0.1 0 0.1 归一化频率序号k 0.2 0.3 0.4 0.5
c) 窗函数 (n) 的长度 L M
X i ( ) DFT [ xi (n) w(n)] xi (n) w(n)e jn
n 0
M 1
S xi ( )


1 | X i ( ) |2 M
1 1 K 从而： S xx ( ) S xi ( ) ，计算流程图可表示如下： U K i 1
Welch法 （ 编 程 ） 求 功 率 谱 密 度 函 数 0 -10 -20 -30 -40 -50 -60 -70 -80 -0.5
归 一 化 功 率 谱 信 号 S-wel /dB
-0.4
-0.3
-0.2
-0.1 0 0.1 归一化频率序号k
wenku.baidu.com
0.2
0.3
0.4
0.5
Welch法 （ 函 数 ） 求 功 率 谱 密 度 函 数 0 -10 -20 -30 -40 -50 -60 -70 -80 -0.5
N 1 j X ( e ) x(n)e jn n 0 2 N 1 X (k ) x(n)e j N kn n 0
S PER ( ) S PER (k )

1 | X (e j ) |2 N 1 | X (k ) |2 N k 0,..., N 1
归 一 化 功 率 谱 信 号 /dB
-30 -40 -50 -60 -70 -80 -90 -100 -0.5
-0.4
-0.3
-0.2
-0.1 0 0.1 归一化频率序号k
0.2
0.3
0.4
0.5
(2) 不同的数据重叠长度对 Welch 谱估计结果有何影响？
不 同 叠 合 长 度 的 Welch谱 估 计 函 数 ， 分 段 长 度 为 70 0 69 35 1