周期图法功率谱估计------窗口效应
FFT和功率谱估计
FFT和功率谱估计1.用Fourier变换求取信号的功率谱---周期图法clf;Fs=1000;N=256;Nfft=256;%数据的长度和FFT所用的数据长度n=0:N-1;t=n/Fs;%采用的时间序列xn=sin(2*pi*50*t)+2*sin(2*pi*120*t)+randn(1,N);Pxx=10*log10(abs(fft(xn,Nfft).^2)/N);%Fourier振幅谱平方的平均值，并转化为dBf=(0:length(Pxx)-1)*Fs/length(Pxx);%给出频率序列subplot(2,1,1),plot(f,Pxx);%绘制功率谱曲线xlabel('频率/Hz');ylabel('功率谱/dB');title('周期图 N=256');grid on;Fs=1000;N=1024;Nfft=1024;%数据的长度和FFT所用的数据长度n=0:N-1;t=n/Fs;%采用的时间序列xn=sin(2*pi*50*t)+2*sin(2*pi*120*t)+randn(1,N);Pxx=10*log10(abs(fft(xn,Nfft).^2)/N);%Fourier振幅谱平方的平均值，并转化为dBf=(0:length(Pxx)-1)*Fs/length(Pxx);%给出频率序列subplot(2,1,2),plot(f,Pxx);%绘制功率谱曲线xlabel('频率/Hz');ylabel('功率谱/dB');title('周期图 N=256');grid on;2.用Fourier变换求取信号的功率谱---分段周期图法%思想：把信号分为重叠或不重叠的小段，对每小段信号序列进行功率谱估计，然后取平均值作为整个序列的功率谱clf;Fs=1000;N=1024;Nsec=256;%数据的长度和FFT所用的数据长度n=0:N-1;t=n/Fs;%采用的时间序列randn('state',0);xn=sin(2*pi*50*t)+2*sin(2*pi*120*t)+randn(1,N);Pxx1=abs(fft(xn(1:256),Nsec).^2)/Nsec; %第一段功率谱Pxx2=abs(fft(xn(257:512),Nsec).^2)/Nsec;%第二段功率谱Pxx3=abs(fft(xn(513:768),Nsec).^2)/Nsec;%第三段功率谱Pxx4=abs(fft(xn(769:1024),Nsec).^2)/Nsec;%第四段功率谱Pxx=10*log10(Pxx1+Pxx2+Pxx3+Pxx4/4);%Fourier振幅谱平方的平均值，并转化为dBf=(0:length(Pxx)-1)*Fs/length(Pxx);%给出频率序列subplot(2,1,1),plot(f(1:Nsec/2),Pxx(1:Nsec/2));%绘制功率谱曲线xlabel('频率/Hz');ylabel('功率谱/dB');title('平均周期图（无重叠） N=4*256');grid on;%运用信号重叠分段估计功率谱Pxx1=abs(fft(xn(1:256),Nsec).^2)/Nsec; %第一段功率谱Pxx2=abs(fft(xn(129:384),Nsec).^2)/Nsec;%第二段功率谱Pxx3=abs(fft(xn(257:512),Nsec).^2)/Nsec;%第三段功率谱Pxx4=abs(fft(xn(385:640),Nsec).^2)/Nsec;%第四段功率谱Pxx5=abs(fft(xn(513:768),Nsec).^2)/Nsec;%第四段功率谱Pxx6=abs(fft(xn(641:896),Nsec).^2)/Nsec;%第四段功率谱Pxx7=abs(fft(xn(769:1024),Nsec).^2)/Nsec;%第四段功率谱Pxx=10*log10(Pxx1+Pxx2+Pxx3+Pxx4+Pxx5+Pxx6+Pxx7/7);%Fourier振幅谱平方的平均值，并转化为dBf=(0:length(Pxx)-1)*Fs/length(Pxx);%给出频率序列subplot(2,1,2),plot(f(1:Nsec/2),Pxx(1:Nsec/2));%绘制功率谱曲线xlabel('频率/Hz');ylabel('功率谱/dB');title('平均周期图（重叠1/2） N=1024');grid on;3.用Fourier变换求取信号的功率谱---welch方法%思想：welch法采用信号重叠分段，加窗函数和FFT算法等计算一个信号序列的自功率谱（PSD）和两个信号序列的互功率谱（CSD），采用MATLAB自%带的函数psdclf;Fs=1000;N=1024;Nfft=256;n=0:N-1;t=n/Fs;window=hanning(256);noverlap=128;dflag='none';randn('state',0);xn=sin(2*pi*50*t)+2*sin(2*pi*120*t)+randn(1,N);Pxx=psd(xn,Nfft,Fs,window,noverlap,dflag);f=(0:Nfft/2)*Fs/Nfft;plot(f,10*log10(Pxx));xlabel('频率/Hz');ylabel('功率谱/dB');title('PSD--Welch方法');grid on;4.功率谱估计----多窗口法（multitaper method ,MTM法）%思想：利用多个正交窗口获得各自独立的近似功率谱估计，综合这些得到一个序列的功率谱估计；相对于普通的周期图有更大的自由度；MTM法采用一个参数：时间带%宽积NW，这个参数用以定义计算功率谱所用窗的数目为2*NW-1，NW越大，时间域分辨率越高而频率分辨率越低，使得功率谱估计的波动减小；随着NW 的增大%，每次估计中谱泄露增多，总功率谱估计的偏差增大clf;Fs=1000;N=1024;Nfft=256;n=0:N-1;t=n/Fs;randn('state',0);xn=sin(2*pi*50*t)+2*sin(2*pi*120*t)+randn(1,N);[Pxx1,f]=pmtm(xn,4,Nfft,Fs); %此处有问题subplot(2,1,1),plot(f,10*log10(Pxx1));xlabel('频率/Hz');ylabel('功率谱/dB');title('多窗口法（MTM）NW=4');grid on;[Pxx,f]=pmtm(xn,2,Nfft,Fs);subplot(2,1,2),plot(f,10*log10(Pxx));xlabel('频率/Hz');ylabel('功率谱/dB');title('多窗口法（MTM）NW=2');grid on;5.功率谱估计----最大熵法（maxmum entmpy method,MEM法）%思想：假定随机序列为平稳高斯过程利用已知的自相关序列rxx(0),rxx(1),rxx(2)...rxx(p)为基础，外推自相关序列rxx(p+1),rxx(p+2)...保证信息熵最大clf;Fs=1000;N=1024;Nfft=256;n=0:N-1;t=n/Fs;window=hanning(256);randn('state',0);xn=sin(2*pi*50*t)+2*sin(2*pi*120*t)+randn(1,N);[Pxx1,f]=pmem(xn,14,Nfft,Fs); %此处有问题subplot(2,1,1),plot(f,10*log10(Pxx1));xlabel('频率/Hz');ylabel('功率谱/dB');title('最大熵法（MEM）Order=14');grid on;%采用Welch方法估计功率谱noverlap=128;dflag='none';subplot(2,1,2)psd(xn,Nfft,Fs,window,noverlap,dflag);xlabel('频率/Hz');ylabel('功率谱/dB');title('Welch方法估计功率谱');grid on;6.功率谱估计----多信号分类法（multiple signal classification,music法）%注：适用于白白噪声中的多正弦波频率估计%思想：将数据自相关矩阵看成是由信号自相关矩阵和噪声自相关矩阵两部分组成，求他们的矩阵特征值向量clf;Fs=1000;N=1024;Nfft=256;n=0:N-1;t=n/Fs;randn('state',0);xn=sin(2*pi*100*t)+2*sin(2*pi*200*t)+randn(1,N);pmusic(xn,[7,1.1],Nfft,Fs,32,16);xlabel('频率/KHz');ylabel('功率谱/dB'); title('Welch方法估计功率谱');grid on;。
经典功率谱和Burg法的功率谱估计
现代信号处理作业实验题目：设信号)()8.0cos(25.0)47.0cos()35.0cos()(321n v n n n n x ++++++=θπθπθπ，其中321,,θθθ是[]ππ,-内的独立随机变量，v(n)是单位高斯白噪声。
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 个观察值x N (n)直接进行傅里叶变换，得到X N (e jw ),然后取其幅值的平方，再除以N ，作为对x （n ）真实功率谱的估计。
2^)(1)(jw e X Nw P N per =， 其中∑-=-=1)()(N n jwn N jwN e n x e X 2）。
BT 法：对于N 个观察值x(0),x(1),。
,x(N-1),令x N (n)=a(n)x(n)。
计算r x （m ）为∑--=-≤+=mN n N Nx N m m n x n xN m r 101),()(1)(,计算其傅里叶变换∑-=--≤=MMm jwm xBT N M e m rm v w P 1 ,)()()(^^，作为观察值的功率谱的估计。
其中v(m)是平滑窗。
3)。
Welch 法：假定观察数据是x(n),n=0,1,2...,N-1，现将其分段，每段长度为M,段与段之间的重叠为M-K,第i 个数据段经加窗后可表示为 1,...,1,0 )()()(-=+=M i iK n x n a n x i M其中K 为一整数，L 为分段数，该数据段的周期图为2)(1)(^w X MU w P i M iper =，其中∑-=-=10)()(M n j w n iM i M e n x w X 。
功率谱估计方法的比较
功率谱估计方法的比较1.周期图法周期图法是最简单直观的功率谱估计方法之一，通过将信号分成多个长为N的区间，计算每个区间内信号的一维傅里叶变换，然后将这些变换结果平方并取平均得到功率谱。
该方法简单快速，但由于其需要使用多个区间的数据进行平均，因此对信号长度有较高的要求，且在信号存在非平稳性时，该方法不适用。
2.自相关法自相关法是一种经典的功率谱估计方法，通过计算信号的自相关函数来估计功率谱。
具体步骤是将信号与其自身的延迟序列进行点乘，并取平均得到自相关函数。
然后对自相关函数进行傅里叶变换，得到功率谱估计值。
该方法计算简单，但精度一般，且在信号长度较长时计算复杂度较高。
3.傅里叶变换法傅里叶变换法是一种经典的功率谱估计方法，通过对信号直接进行傅里叶变换得到功率谱。
该方法计算简单，精确度高，但对信号的长度存在要求，较长的信号长度能提供更高的分辨率。
此外，傅里叶变换法只适用于周期性信号。
4.平均周期图法平均周期图法是一种对周期图法的改进。
它将信号分为多段，并对每一段进行周期图计算，然后将计算结果平均得到平均周期图。
与周期图法相比，平均周期图法可以降低误差，提高估计精度。
然而，该方法仍然对信号长度有一定要求，并且计算复杂度较高。
5.移动平均法移动平均法是一种基于滑动窗口的功率谱估计方法，其基本思想是通过对信号进行多次滑动窗口处理，将窗口内信号的傅里叶变换结果平方并取平均得到功率谱估计值。
该方法在计算复杂度上较低，适用于非平稳信号的功率谱估计。
但是，由于窗口大小的选择存在权衡，需要根据实际情况进行合理设置。
总结起来，各种功率谱估计方法各有优劣。
周期图法和自相关法计算简单，但方法的精度较低，受信号长度限制且无法处理非平稳信号。
傅里叶变换法具有较高的计算精度，但对信号的长度和周期性要求较高。
平均周期图法和移动平均法对周期图法进行了改进，在精度上有所提高，但计算复杂度较高。
因此，在实际应用中，需要根据具体的信号特点和处理要求选取合适的功率谱估计方法。
功率谱估计的方法
功率谱估计的方法
功率谱估计是信号处理中常用的一种方法，用于分析信号在频域内的特点，通常可以分为以下几种方法：
一、经典方法
1.傅里叶变换法：将时域信号通过傅里叶变换变换到频域，然后计算功率谱密度。
2.自相关法：通过自相关函数反映信号的统计平稳性，然后通过傅里叶变换计算功率谱密度。
3.周期图法：将信号分解为若干个周期波形，然后对每个周期波形进行傅里叶变换计算周期功率谱，最后汇总得到整个信号的功率谱。
二、非经典方法
1. 时-频分析法：如短时傅里叶变换（STFT）、小波变换等，将信号分解为时域和频域两个维度的分量，从而可以分析信号在时间和频率上的变化。
2. 基于协方差矩阵的特征值分解法：通过建立协方差矩阵，在张成空
间中求解特征向量，从而达到计算信号功率谱的目的。
3. 基于频率锁定法：如MUSIC法、ESPRIT法等，是一种利用特定信号空间中的特定模式进行处理的方法。
以上方法各有特点，根据实际需求选择不同的方法可以得到相应的功率谱估计结果。
生医信号处理课件6经典功率谱
（2）方差 Var[S (e j )] Var[S (e j )] S 2 (e j ) M per x 频率分辨率及旁瓣泄漏
Re s[SM (e j )] ()3dB
二、平均法（Bartlett法）
Bartlett提出，将 x(n)分为长L，互不重叠的k段子序列，N=kL,
3、讨论
m M 2 w ( m) M
M↑偏差↓
M↑方差↑ 一般M=N/5
五、谱估计技术的性能指标
1、变异性γ （归一化方差）
2、品质因数μ
4 2 x
j
j
结论：非一致估计
例
三、周期图的随机起伏
Cov[S per (e
j1
), S per (e
j2
sin[(k l ) ] 2 )] [ ] (k l ) N sin[ ] N
4 x
一、数据加窗（修正周期图）
数据窗
6.4.3 功率谱估计的改进
E[ S per (e )]
j
j
1 E[ S BT (e )] S x ( e j ) W ( e j ) 2
j
N ,
W (e j ) ()
是渐近无偏
2、方差
if 1 M N
j
1 2 Var[ S BT (e )] S x (e j ) W 2 (e j ( u ) )du 2N 1 j S x (e ) N
Bartlett法和Welch法分别对周期图和修正周期图进行平均， 从而达到减少方差的目的。 Blackman-Tukey法为了减少周
期图的方差，对自相关序列的估计进行加窗处理，从而减
少自相关序列的估计中那些不可靠的估计值对周期图的贡 献。
周期图法功率谱估计------窗口效应
-35
50
300
窗函数对估计结果的影响主要取决于其主瓣 和旁瓣特性，所以窗函数评价指标有 。
主瓣宽度
窗函 数评 价指 标
旁瓣大小
旁瓣衰减斜率
谱估计法比较
周期图法功率谱估计其特点是离散性大，曲 线粗糙，方差较大，但分辨率较高。 窗函数周期图法功率谱估计的收敛性较好, 曲 线平滑, 估计的结果方差较小, 但是功率谱主 瓣较宽, 分辨率低。这是由于对随机序列的分 段处理引起了长度有限所带来的Gibbs 现象 而造成的。 窗口效应的谱估值比较平滑, 但是分辨率较差 。其原因是给每一段序列用适当的窗口函数 加权后, 在得到平滑的估计结果的同时, 使功 率谱的主瓣变宽, 因此分辨率有所下降。
MALAB中可以采用下面的命令来生成：
数据分段 窗处理 各段功率谱
平均功率谱
Welch法谱估计流程图
Welch 法优点, 一是选择适当的窗函数w(n), 并在周期 图计算前直接加进去, 加窗的优点是无论什么样的窗函 数均可使谱估计非负。二是在分段时, 可使各段之间有 重叠,这样会使方差减小。
加窗后的窗口效应仿真
-5
-10
-15
-20
-25
-30
-35
0
50
100
150
200
250
300
海明窗处理谱估计 的旁瓣部分衰减较 大，方差较小，噪 声水平较低，性能 良好，改善了由矩 形窗处理的谱估计 所产生的较大谱失 真
加布莱克曼窗后的窗口效应
加 blackman 窗 0
-5
-10
-15
-20
-25
-30
布莱克曼窗处理谱 估计的旁瓣部分衰 减较大，方差较小， 噪声水平较低，性 能良好，改善了由 矩形窗处理的谱估 计所产生的较大谱 失真
改进周期图法功率谱估计中的窗函数仿真分析
仿真 ， 就其中的频率分辨率 、 采样信号的数据长度与窗函数的特性等 因素 之间的关系 ， 进行 了详细的讨论 ， 重点分析 了窗函 数对谱估计 的影 响，指出了它们 的优缺点， 出了不 同数据类 型与估计方 式下 的功 率谱标 准方差 的均值 。最后 ， 得 根据 实验 结果并结合 窗函数 的评价指标 ， 出了选取合适的窗函数应注意的一般原则 ， 指 具有 一定 的实际指导意义 。 关键词 ： 周期 图； 谱估计 ； 函数 ； 窗 韦尔奇法
ｂ ｓ ｄ ｏ ｅｃ ．Ｍｅ ｎ ｈ ｌ ｔｐ ｉ ｒ ｙ ａ ａｙ ｅ ｈ ｎ ｌｅ ｃ ｆｗ ｎ ｏ ｆ ｎ ｔ ｎ ｎ ｏｎ ｓｏ ｔ ｈ ｉ ｖｒ ｅ ｄ ａ ｅ ｎ Ｗ ｌｈ ａｗ ｉ ｅ，ｉ ｒｍａ ｉ ｎ ｌｚ ｓｔ ｅ ｉｆｕ ｎ ｅ ｏ ｄ ｗ ｃｉ ｓ ａ ｄ ｐ ｉｔ ｕ ｅｒ ｉ ｕ ｓ ａ ｌ ｉ ｕ ｏ ｔ ｔ ｎ
（ ．重庆大学通信与测控研究所 ， １ 重庆 ４ ０ ４ ；．重庆通信学院 ， ０ ０４ ２ 重庆 ４ ０ ３ ） ０ ０ ５ 摘要 ： 图法是功率谱估计 的一种基本方法 ， 周期 但该方法不满足一致 性估计 条件 ， 谱估计 的分辨力和方差都很难满 足实际应 用需要。因此采用其改进方法 （ Ｗ ｌ 即 ｅ ｈ法） ｃ 估计信号功率谱 ， 并结合该方法对四种典型窗函数的谱估计性能进行 了具体的
Ｙｕ ｎ—ｆ ｎ ， Ｘｕ ｅ ｇ ＭＡ Ｄａ—ｗｅ ， ＥＩＬｉ ｉ Ｗ ｎ
窗口效应
一种改进的窗处理是韦尔奇（Welch）法，因为窗在两 边渐变为零，所以这种方法降低了由于重叠导致的段间 统计相依的效应，而且，一个合适的非矩形窗可以减小 “旁瓣效应”，即“谱泄漏”。 经常采用的窗函数有矩形窗、汉宁窗(Hanning)、海明窗 (Hamming)、凯瑟窗(Kaiser)和切比雪夫窗(Chebwin)等， 其中矩形窗、汉宁窗以及海明窗都是广义余弦窗的特殊 形式，它们可以看作是频率为0、 ／(N一1)和 ／(N一1)的余弦曲线的线性组合，这里N为窗的长度，
巴特利特平均周期图法的实现过程就是把数据分成互 不重叠的L 段，每段有M个样本，即
其中 令
, 是长度为M的矩形窗。 ,则第i段的周期图可有下面公式求得：
由上式可求得整个序列的谱估计：
巴特利特平均周期图法的仿真
在巴特利特平均周期图 法中，对于固定长度的 N，分段K数目愈大， 方差愈小，但M愈小， 会使偏倚增大，谱分辨 率变得愈差。 因此在应用巴特利特平 均周期图法时，需要在 偏倚和估计方差之间进 行权衡，根据实际要求 选择适当的M和N。
在MALAB中可以采用下面的命令来生成：
加窗后的窗口效应
加矩形窗 0
-5
-10
-15
-20
-25
-30
-35
0
50
100
150
200
250
300
矩形窗处理的谱估 计的主瓣较窄，分 辨率最好，同时其 主瓣附近的衰减比 其他两个窗函数小， 旁瓣部分泄漏比较 严重，所以其方差 特性最差，噪声水 平较高。
例如测量震动物体的自震频率时，则可以选用主瓣宽度比 较窄的Rectangular窗；当要求谱泄漏小时，则应选用旁 瓣幅度较小的窗。
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
一种改进的窗处理是韦尔奇（Welch）法，因为窗在两 边渐变为零，所以这种方法降低了由于重叠导致的段间 统计相依的效应，而且，一个合适的非矩形窗可以减小 “旁瓣效应”，即“谱泄漏”。 经常采用的窗函数有矩形窗、汉宁窗(Hanning)、海明窗 (Hamming)、凯瑟窗(Kaiser)和切比雪夫窗(Chebwin) 等，其中矩形窗、汉宁窗以及海明窗都是广义余弦窗的 特殊形式，它们可以看作是频率为0、 2 ( N 1) 和 4 ( N 1) 的余弦曲线的线性组合，这里N为窗的长度，在
jw
注：1.做傅立叶变换时，为方便计算机处理，通常将数据补足2^n点 2.FFT程序的序号K是归一频率序号，转换为实际频率为k*Fs/N
周期图法仿真
Fs=600; nfft=512; n=0:1/Fs:1; xn=cos(2*pi*40*n)+3*cos(2*pi*90*n)+randn(size(n)); X=fft(x,nfft); Pxx=abs(X).^2/length(n); %求解PSD t=0:round(nfft/2-1); f=t*Fs/nfft; P=10*log10(Pxx(t+1)); %纵坐标的单位为dB plot(f,P)
Pxx ( w)
M 1 m ( M 1) jwm ( m ) w ( m ) e xx
式中窗函数序列w(m) 的长度为2M-1的有限宽 窗序列，为了保证 Pxx (w) 是一个实、偶、非负 函数， w(m)应是一个偶序列，并且满足条件
W (e ) 0, w
加矩形窗后的窗口效应
加矩形窗 0
-5
-10
-15
-20
-25
-30
-35
0
50
100
150
200
250
300
矩形窗处理的谱估 计的主瓣较窄，分 辨率最好，同时其 主瓣附近的衰减比 其他两个窗函数小， 旁瓣部分泄漏比较 严重，所以其方差 特性最差，噪声水 平较高。
加海明窗后的窗口效应
加海明窗 0
-35
0
50
100
150
200
250
300
窗函数对估计结果的影响主要取决于其主瓣 和旁瓣特性，所以窗函数评价指标有 。
主瓣宽度
窗函 数评 价指 标
周期图法功率谱估计------窗口效应
主要内容
概述 周期图法 窗函数法
谱估计法比较
总结
数字信号处理的两个主要分支：
数字滤波 频谱分析
对随机信号的频谱分析——功率谱估计
对确定信号，可以用傅立叶变换；而随机 信号无始无终具无限能量，不满足傅立叶变 换绝对可积的条件.
概述
功率谱：随机信号的功率谱反映的是随机 信号的频率成分及各成分的相对强弱。 功率谱估计：基于有限的数据寻找信号、 随机过程或系统的频率成分。
两种基本算法： 周期图法：把功率谱和信号幅频特性的平 方结合起来。 自相关法： 根据维纳-辛钦定理，先估计相 关函数，再经傅立叶变换得功率谱估计。
周期图法
X e jw
^
n jnw jnw x e x e n n n 0 N 1
1 jw 2 S x e X e N
-5
-10
-15
-20
-25
-30
-35
0
50
100
150
200
250
300
海明窗处理谱估计 的旁瓣部分衰减较 大，方差较小，噪 声水平较低，性能 良好，改善了由矩 形窗处理的谱估计 所产生的较大谱失 真
加布莱克曼窗后的窗口效应
加 blackman 窗 0
-5
-10
-15
-20
-25
-30
布莱克曼窗处理谱 估计的旁瓣部分衰 减较大，方差较小， 噪声水平较低，性 能良好，改善了由 矩形窗处理的谱估 计所产生的较大谱 失真
•plot_Pxx=10*log10(Pxx); •plot_Pxx1=10*log10(Pxx1); •plot_Pxx2=10*log10(Pxx2); •figure(1) •plot(f,plot_Pxx); •title('加矩形窗'); •figure(2) •plot(f,plot_Pxx1); •title('加海明窗'); •figure(3) •plot(f,plot_Pxx2); •title('加blackman 窗');
Matlab 代码: Fs=600; n=0:1/Fs:1; xn=cos(2*pi*40*n)+3*cos(2*pi*90*n)+randn(size (n)); nfft=512; window=boxcar(100);%矩形窗 window1=hamming(100);%海明窗 window2=blackman(100);%blackman 窗 noverlap=20; %数据无重叠 range='onesided'; %频率间隔为[0 Fs/2], 计算一半 的频率 [Pxx,f]=pwelch(xn,window,noverlap,nfft,Fs,range ); [Pxx1,f]=pwelch(xn,window1,noverlap,nfft,Fs,ran ge); [Pxx2,f]=pwelch(xn,window2,noverlap,nfft,Fs,ran ge);
MALAB中可以采用下面的命令来生成：
数据分段 窗处理 各段功率谱
平均功率谱
Welch法谱估计流程图
Welch 法优点, 一是选择适当的窗函数w(n), 并在周期 图计算前直接加进去, 加窗的优点是无论什么样的窗函 数均可使谱估计非负。二是在分段时, 可使各段之间有 重叠,这样会使方差减小。
加窗后的窗口效应仿真
Байду номын сангаас 周期图法仿真
30 20
数据长度N太大， 谱线起伏加剧
10
0
-10
-20
数据长度N太小， 谱的分辨率不 好
0 50 100 150 200 250 300
-30
功率谱估计的改进——窗函数法
加窗处理的宗旨是减小频谱泄漏，但却降低了 频谱分辨率。一般说来，减小泄漏和提高分辨 率是互相矛盾的，但是如何使两者都取得最好 的效果，就必须对窗函数的特点进行分析，然 后根据实际要求来选择。其公式如下：