归 一 化 功 率 谱 信 号 /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
-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 可以看出，自相关法得到的结 果与周期图法相似，同样是方差较大，信号振荡大。
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
0.2
0.3
0.4
0.5
Welch法 （ 函 数 ） 求 功 率 谱 密 度 函 数 0 -10 -20 -30 -40 -50 -60 -70 -80 -0.5
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 为得到 的频率点。 分别采用编程和函数两种方法求解所有数据的功率谱， 并且进行归一化后的 结果如下图所示：
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
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
加窗平滑法求功率谱 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 段
功率谱估计分析及比较
1 实验目的
(1) 掌握Welch算法的概念、应用及特点； (2) 了解谱估计在信号分析中的作用； (3) 能够利用 Welch 法对信号作谱估计，对信号的特点加以分析。
2 实验内容
(1) 读入实验数据。 (2) 编写一利用Welch法作估计的算法程序。 (3) 将计算结果表示成图形的形式，给出信号谱的分布情况图。
3 算法讨论与分析
随机信号的功率密度函数的定义如下：
S xx ( j) lim E[| X T ( j) |2 ]
T
其物理意义表示为信号在单位频带内的平均功率。 经典的谱估计主要包括以 下几种方式： (1) 周期图法 周期图法是直接建立在功率谱的定义式上的，也称之为直接法。原理计算如 下：
对于直接法的功率谱估计，当数据长度 N 太大时，谱曲线起伏加剧，若 N 太小，谱的分辨率又不好，Welch 谱估计法从两个方面对直接法进行改进： 一是分段，且各段之间可以存在区间重叠，这样会使方差减小； 二是加窗， 选择适当的窗函数加在 x(n) 的分段数据上， 这样可使谱估计变得 更加平滑，但分辨率减小。 从理论分析上来说， Welch 谱估计法是渐进无偏估计， 并且是渐进一致估计。
(3) 不同的窗函数对 Welch 谱估计结果有何影响？
归 一 化 功 率 谱 信 号 S /dB
不 同 窗 函 数 下 的 Welch谱 估 计 函 数 ， 叠 合 长 度 为 50% ， 每 段 长 度 为 100 0 矩形窗 -10 三角窗 布莱克曼窗 汉明窗 -20 汉宁窗 -30 -40 -50 -60 -70 -80 -90 -0.5
归 一 化 功 率 谱 信 号 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)
自相关法求功率谱密度函数 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 -20 -40 0 -20 -40 直接法（函数）求功率谱密度函数