经典功率谱设计

DSP 实验报告

实验题目: 实验 功率谱估计

一、实验要求：

（1）理解功率谱估计的基本概念；

（2）掌握经典功率谱估计方法——直接法和间接法； （3）掌握改进的经典功率谱估计方法，例如Welch 法。

二、实验内容与原理：

功率谱估计就是基于有限的数据寻找信号、随机过程或系统的频率成分。它表示随机信号频率域的统计特性。随机信号是无始无终具有无限能量的，所以其傅立叶变换并不存在，因为它不满足绝对可积的条件。因此需要研究其在频率域上的功率分布情况，即功率谱密度或功率谱。

根据实验要求，完成该实验首先要正确的生成被估计信号x (n ) 。x (n ) 数据的长度和FFT 所用的数据长度都设为1024。

1. 周期图法：

直接法，即周期图法，是由傅立叶变换得到的：将随机信号x (n ) 的N 点样本值x N (n ) 看作能量有限信号，取其傅立叶变换，得到X N (e jw ) ；然后再取其幅值的平方，并除以N 作为x (n ) 的真实功率谱P (e jw ) 的估计，即

P (e jw ) =

1

|X N (w ) |2 N

实验中，将随机信号x(n)的N 点样本值看作xN(n)能量有限信号，取其傅立叶变换，得到X ，然后再取其幅值的平方，并除以N 作为x(n)的真实功率谱P 的估计。

2. 间接法：

间接法，又称为自相关法或BT 法，是由随机信号N 个观察值

x (0),..., x (N -1) ，估计出自相关函数R N (m ) ，然后再求R N (m ) 的傅立叶变换作为

功率谱的估计：

S (e ) =

jw

m =-M

∑R

M

N

(m ) e -jwm |M |≤N -1

即如下计算：

1N -1

ˆx (m ) =∑x N (n ) x N (n +m ) r

N n =0

1=N

N -1-m

n =0

x N (n ) x N (n +m )

|m |≤N -1

ˆ(ω) =P BT

m =-M

∑v (m ) r ˆ(m ) e

x

M

-j ωm

M ≤N -1

实验中由随机信号N 个观察值估计出自相关函数R(m),然后再求R(m)傅立叶变换作为功率谱的估计PBT 。

直接法和间接法的方差性能很差，而且当数据长度太大时，谱曲线起伏加剧；若数据长度太小，则谱的分辨率又不好，所以需要改进。改进的直接谱估计方法由Bartlett 法和Welch 法。

3.BARTLETT 算法

Bartlett 法将采样数据x N (n ) 分成L 段，每段的长度都是M ，即N=LM，对

i

每段数据加矩形窗，再计算其各自的功率谱P PER (w ) ，把P P i E R (w ) 对应相加，再取

_

平均，得到平均周期图P PER (w ) 。即如下过程：

首先将观测数据分为L 段，每段长M ，分段后的数据可以用式(1)表示：

i

x M (n ) =x (n +iM )

i =0,1,..., L -1 ；n =0,1,..., M -1；LM ≤N (1)

其中

i X (ω) =∑x M (n ) e -j ωn i M

n =0M -1

平均的周期图为：

1L -1ˆi

P per (ω) ＝∑P per (ω)

L i =0

_

BARTLETT 法是周期图算法的一种改进。由概率论的知识知道，如果

X 1, X 2 X N 是N 个不相关的随机变量，每个随机变量的期望值为μ，方差为

σ2，那么将这N 个随机变量求平均，它的期望仍为μ，方差变为σ。BARTLETT 法即是受此启发，将观测数据分段，先求每段数据的周期图，再求平均的周期图，当分段较多时，估计出的功率谱较平滑，频率分辨率较差；当分段较少时，估计出的功率谱起伏较大，频率分辨率较好。

4.WELCH 算法

Welch 法是对Bartlett 法的改进：一，在对x N (n ) 分段时，可允许每段数据有部分重叠；二，每段数据窗口可以不是矩形窗口，例如使用汉宁窗或哈明窗，记为d (n ) 。然后按

Bartlett

法求每一段的功率谱，记为

M -1n =0

2

P

i PER

1M -1i 1|∑x N (n ) d (n ) e -jwn |2，其中U =(w ) =

MU n =0M

N -M /2

）

M /2

∑d

2

(n ) 。平均后的功率谱为：

（段数L =

1L i

P P E (P P E (∑R w ) =R w )

L i =1

加窗的优点是使得无论对于什么样的窗函数均可以谱估计为非负值；二是在分段时，各段之间有重叠，这样会使方差减小。

三、实验目的：

分别用四种不同的发放进行功率谱估计，并对比结果。

四、实验程序：

clear;

%数据的长度和FFT 所用的数据长度 nfft=1024;N=1024; %每段长度 Ns=256;

%产生含有噪声的序列xn n=[0:N-1];

w1=2*pi*0.02; w2=2*pi*0.28; wn=randn(1,N);

xn=sin(w1*n)+2*cos(w2*n)+wn; %直接法求功率谱

%将随机信号x(n)的N 点样本值看作xN(n)能量有限信号，取其傅立叶变换，得到X ；然后再取其幅值的平方，并除以N 作为x(n)

%的真实功率谱P 的估计。 %计算序列的DFT XN=fft(xn,nfft);

%对序列取绝对值后平方 PER=abs(XN).^2/N; %并转化为dB PERdb=10*log10(PER); %给出频率序列

f=(0:length(PERdb)-1)/length(PERdb); %绘制功率谱图形 figure(1); plot(f,PERdb); xlabel('频率/Hz'); ylabel('功率谱/dB'); title('直接法 N=1024'); grid;

%间接法求功率谱

%又称为自相关法或BT 法，是由随机信号N 个观察值估计出自相关函数R(m),然后再求R(m)傅立叶变换作为功率谱的估计PBT

%计算序列的自相关函数Rm Rm=xcorr(xn,'unbiased');

%计算自相关函数Rm 的DTFT PBT=fft(Rm,nfft); %把PBT 转化为dB PBTdb=10*log10(abs(PBT));

Fbt=(0:length(PBTdb)-1)/length(PBTdb); %绘制功率谱图形 figure(2); plot(Fbt,PBTdb); xlabel('频率/Hz'); ylabel('功率谱/dB'); title('间接法 N=1024'); grid;

%bartlett法求功率谱

%Bartlett平均周期图的方法是将N 点的有限长序列x(n)分段，对各段用周期图法求解功率后再平均。

%加矩形窗 window=ones(1,Ns); normvalue=norm(window);

PBAR1=abs(fft(window.*xn(1:256),Ns).^2)/normvalue^2;%第一段功率谱 PBAR2=abs(fft(window.*xn(257:512),Ns).^2)/normvalue^2;%第二段功率谱 PBAR3=abs(fft(window.*xn(513:768),Ns).^2)/normvalue^2;%第三段功率谱 PBAR4=abs(fft(window.*xn(769:1024),Ns).^2)/normvalue^2;%第四段功率谱 %求Fourier 振幅谱的平均值，并转化为dB

PBAR=10*log10((PBAR1+PBAR2+PBAR3+PBAR4)/4); %给出频率序列

Fpbar=(0:length(PBAR)-1)/length(PBAR); %绘制功率谱曲线 figure(3);

plot(Fpbar,PBAR);

xlabel('频率/Hz');ylabel('功率谱/dB'); title('bartlett法 4*256'); grid;

%welch法求功率谱

%Welch法对Bartlett 法进行了两方面的修正，一是选择适当的窗函数w(n)，并再周期图计算前直接加进去，加窗的优点是无论什么样的窗函数均可使谱估计非负。

%二是在分段时，可使各段之间有重叠，这样会使方差减小。 %加汉宁窗

hanningwindow=hanning(Ns)'; normvalue=norm(window);

PWel1=abs(fft(hanningwindow.*xn(1:256),Ns).^2)/normvalue^2;%第一段功率谱

PWel2=abs(fft(hanningwindow.*xn(129:384),Ns).^2)/normvalue^2;%第二段功率谱

PWel3=abs(fft(hanningwindow.*xn(257:512),Ns).^2)/normvalue^2;%第三段功率谱

PWel4=abs(fft(hanningwindow.*xn(385:640),Ns).^2)/normvalue^2;%第四段功率谱

PWel5=abs(fft(hanningwindow.*xn(513:768),Ns).^2)/normvalue^2;%第五段功率谱

PWel6=abs(fft(hanningwindow.*xn(641:896),Ns).^2)/normvalue^2;%第六段功率谱

PWel7=abs(fft(hanningwindow.*xn(769:1024),Ns).^2)/normvalue^2;%第七段功率谱

%求Fourier 振幅谱的平均值，并转化为dB

PWel=10*log10((PWel1+PWel2+PWel3+PWel4+PWel5+PWel6+PWel7)/7);

%给出频率序列

FWel=(0:length(PWel)-1)/length(PWel); %绘制功率谱曲线 figure(4); plot(FWel,PWel);

xlabel('频率/Hz');ylabel('功率谱/dB'); title('welch法 7*256'); grid;

五、实验结果：

DSP 实验报告

实验题目: 实验 功率谱估计

一、实验要求：

（1）理解功率谱估计的基本概念；

（2）掌握经典功率谱估计方法——直接法和间接法； （3）掌握改进的经典功率谱估计方法，例如Welch 法。

二、实验内容与原理：

功率谱估计就是基于有限的数据寻找信号、随机过程或系统的频率成分。它表示随机信号频率域的统计特性。随机信号是无始无终具有无限能量的，所以其傅立叶变换并不存在，因为它不满足绝对可积的条件。因此需要研究其在频率域上的功率分布情况，即功率谱密度或功率谱。

根据实验要求，完成该实验首先要正确的生成被估计信号x (n ) 。x (n ) 数据的长度和FFT 所用的数据长度都设为1024。

1. 周期图法：

直接法，即周期图法，是由傅立叶变换得到的：将随机信号x (n ) 的N 点样本值x N (n ) 看作能量有限信号，取其傅立叶变换，得到X N (e jw ) ；然后再取其幅值的平方，并除以N 作为x (n ) 的真实功率谱P (e jw ) 的估计，即

P (e jw ) =

1

|X N (w ) |2 N

实验中，将随机信号x(n)的N 点样本值看作xN(n)能量有限信号，取其傅立叶变换，得到X ，然后再取其幅值的平方，并除以N 作为x(n)的真实功率谱P 的估计。

2. 间接法：

间接法，又称为自相关法或BT 法，是由随机信号N 个观察值

x (0),..., x (N -1) ，估计出自相关函数R N (m ) ，然后再求R N (m ) 的傅立叶变换作为

功率谱的估计：

S (e ) =

jw

m =-M

∑R

M

N

(m ) e -jwm |M |≤N -1

即如下计算：

1N -1

ˆx (m ) =∑x N (n ) x N (n +m ) r

N n =0

1=N

N -1-m

n =0

x N (n ) x N (n +m )

|m |≤N -1

ˆ(ω) =P BT

m =-M

∑v (m ) r ˆ(m ) e

x

M

-j ωm

M ≤N -1

实验中由随机信号N 个观察值估计出自相关函数R(m),然后再求R(m)傅立叶变换作为功率谱的估计PBT 。

直接法和间接法的方差性能很差，而且当数据长度太大时，谱曲线起伏加剧；若数据长度太小，则谱的分辨率又不好，所以需要改进。改进的直接谱估计方法由Bartlett 法和Welch 法。

3.BARTLETT 算法

Bartlett 法将采样数据x N (n ) 分成L 段，每段的长度都是M ，即N=LM，对

i

每段数据加矩形窗，再计算其各自的功率谱P PER (w ) ，把P P i E R (w ) 对应相加，再取

_

平均，得到平均周期图P PER (w ) 。即如下过程：

首先将观测数据分为L 段，每段长M ，分段后的数据可以用式(1)表示：

i

x M (n ) =x (n +iM )

i =0,1,..., L -1 ；n =0,1,..., M -1；LM ≤N (1)

其中

i X (ω) =∑x M (n ) e -j ωn i M

n =0M -1

平均的周期图为：

1L -1ˆi

P per (ω) ＝∑P per (ω)

L i =0

_

BARTLETT 法是周期图算法的一种改进。由概率论的知识知道，如果

X 1, X 2 X N 是N 个不相关的随机变量，每个随机变量的期望值为μ，方差为

σ2，那么将这N 个随机变量求平均，它的期望仍为μ，方差变为σ。BARTLETT 法即是受此启发，将观测数据分段，先求每段数据的周期图，再求平均的周期图，当分段较多时，估计出的功率谱较平滑，频率分辨率较差；当分段较少时，估计出的功率谱起伏较大，频率分辨率较好。

4.WELCH 算法

Welch 法是对Bartlett 法的改进：一，在对x N (n ) 分段时，可允许每段数据有部分重叠；二，每段数据窗口可以不是矩形窗口，例如使用汉宁窗或哈明窗，记为d (n ) 。然后按

Bartlett

法求每一段的功率谱，记为

M -1n =0

2

P

i PER

1M -1i 1|∑x N (n ) d (n ) e -jwn |2，其中U =(w ) =

MU n =0M

N -M /2

）

M /2

∑d

2

(n ) 。平均后的功率谱为：

（段数L =

1L i

P P E (P P E (∑R w ) =R w )

L i =1

加窗的优点是使得无论对于什么样的窗函数均可以谱估计为非负值；二是在分段时，各段之间有重叠，这样会使方差减小。

三、实验目的：

分别用四种不同的发放进行功率谱估计，并对比结果。

四、实验程序：

clear;

%数据的长度和FFT 所用的数据长度 nfft=1024;N=1024; %每段长度 Ns=256;

%产生含有噪声的序列xn n=[0:N-1];

w1=2*pi*0.02; w2=2*pi*0.28; wn=randn(1,N);

xn=sin(w1*n)+2*cos(w2*n)+wn; %直接法求功率谱

%将随机信号x(n)的N 点样本值看作xN(n)能量有限信号，取其傅立叶变换，得到X ；然后再取其幅值的平方，并除以N 作为x(n)

%的真实功率谱P 的估计。 %计算序列的DFT XN=fft(xn,nfft);

%对序列取绝对值后平方 PER=abs(XN).^2/N; %并转化为dB PERdb=10*log10(PER); %给出频率序列

f=(0:length(PERdb)-1)/length(PERdb); %绘制功率谱图形 figure(1); plot(f,PERdb); xlabel('频率/Hz'); ylabel('功率谱/dB'); title('直接法 N=1024'); grid;

%间接法求功率谱

%又称为自相关法或BT 法，是由随机信号N 个观察值估计出自相关函数R(m),然后再求R(m)傅立叶变换作为功率谱的估计PBT

%计算序列的自相关函数Rm Rm=xcorr(xn,'unbiased');

%计算自相关函数Rm 的DTFT PBT=fft(Rm,nfft); %把PBT 转化为dB PBTdb=10*log10(abs(PBT));

Fbt=(0:length(PBTdb)-1)/length(PBTdb); %绘制功率谱图形 figure(2); plot(Fbt,PBTdb); xlabel('频率/Hz'); ylabel('功率谱/dB'); title('间接法 N=1024'); grid;

%bartlett法求功率谱

%Bartlett平均周期图的方法是将N 点的有限长序列x(n)分段，对各段用周期图法求解功率后再平均。

%加矩形窗 window=ones(1,Ns); normvalue=norm(window);

PBAR1=abs(fft(window.*xn(1:256),Ns).^2)/normvalue^2;%第一段功率谱 PBAR2=abs(fft(window.*xn(257:512),Ns).^2)/normvalue^2;%第二段功率谱 PBAR3=abs(fft(window.*xn(513:768),Ns).^2)/normvalue^2;%第三段功率谱 PBAR4=abs(fft(window.*xn(769:1024),Ns).^2)/normvalue^2;%第四段功率谱 %求Fourier 振幅谱的平均值，并转化为dB

PBAR=10*log10((PBAR1+PBAR2+PBAR3+PBAR4)/4); %给出频率序列

Fpbar=(0:length(PBAR)-1)/length(PBAR); %绘制功率谱曲线 figure(3);

plot(Fpbar,PBAR);

xlabel('频率/Hz');ylabel('功率谱/dB'); title('bartlett法 4*256'); grid;

%welch法求功率谱

%Welch法对Bartlett 法进行了两方面的修正，一是选择适当的窗函数w(n)，并再周期图计算前直接加进去，加窗的优点是无论什么样的窗函数均可使谱估计非负。

%二是在分段时，可使各段之间有重叠，这样会使方差减小。 %加汉宁窗

hanningwindow=hanning(Ns)'; normvalue=norm(window);

PWel1=abs(fft(hanningwindow.*xn(1:256),Ns).^2)/normvalue^2;%第一段功率谱

PWel2=abs(fft(hanningwindow.*xn(129:384),Ns).^2)/normvalue^2;%第二段功率谱

PWel3=abs(fft(hanningwindow.*xn(257:512),Ns).^2)/normvalue^2;%第三段功率谱

PWel4=abs(fft(hanningwindow.*xn(385:640),Ns).^2)/normvalue^2;%第四段功率谱

PWel5=abs(fft(hanningwindow.*xn(513:768),Ns).^2)/normvalue^2;%第五段功率谱

PWel6=abs(fft(hanningwindow.*xn(641:896),Ns).^2)/normvalue^2;%第六段功率谱

PWel7=abs(fft(hanningwindow.*xn(769:1024),Ns).^2)/normvalue^2;%第七段功率谱

%求Fourier 振幅谱的平均值，并转化为dB

PWel=10*log10((PWel1+PWel2+PWel3+PWel4+PWel5+PWel6+PWel7)/7);

%给出频率序列

FWel=(0:length(PWel)-1)/length(PWel); %绘制功率谱曲线 figure(4); plot(FWel,PWel);

xlabel('频率/Hz');ylabel('功率谱/dB'); title('welch法 7*256'); grid;

五、实验结果：


相关文章

  • 小松鼠燃气壁挂炉

    • 一.产品名称:M2 蓝钻分段燃烧系列 (先进分段燃烧技术,舒适洗浴享受 ) 品牌:小松鼠燃气壁挂炉 气种:天然气 液化气 人工煤气 型号:采暖/ 生活热水两用,18 / 24 KW 输出功率 产品特点: • 时尚弧线外观设计,尽显现代简约之 ...查看


  • 随机过程课程设计模版(新)

    • <随机过程> 课程设计(论文) 题 目: 平 稳 时 间 序 列 的 MA(q) 模 型 的 预 报 学 院: 理 学 院 专 业: 数 学 与 应 用 数 学 班 级: 数 学 09-1 班 学 生 姓 名: 学 生 学 号: ...查看


  • 数字信号处理试卷

    • 一. 填空(2分/空,共30分) 1. 对一个1Hz 的正弦波信号进行10Hz 抽样.请问该信号的连续角频率Ω是[2πrad/s],圆频率ω是 [0.2πrad ]. 2. 假定信号的功率为P S , 噪声功率为P U , 若信噪比SNR= ...查看


  • 中国工程师标配十款经典RF射频器件

    • 中国工程师标配十款经典RF射频器件 RF射频器件对非射频领域的工程师来说并没有那么悬乎,说白了就是一种可发生高频交流变化电磁波的器件,目前的应用非常广,如智能手机.GPS.手持设备等等. 下面这篇文章肯定不会让你失望,让你更清楚的了解RF射 ...查看


  • 现代信号处理新方法试题

    • 现代信号处理新方法试题 一. 填空题 1.平稳随机信号是指:. 判断随机信号是否广义平稳的三个条是: . 高斯白噪声信号是指: . 信号的遍历性是指: . 广义遍历信号x(n)的时间均值的定义为: 其时间自相关函数的定义为: 件 2.离散随 ...查看


  • 随机信号分析理论的应用综述

    • 随机信号分析理论的应用综述 (结课论文) 学院: 系别:电子信息工程 班级: 姓名: 学号: 指导老师: 目录 第一章 概述 1.1 随机信号分析的研究背景 1.2 随机信号分析的主要研究问题 第二章 随机信号分析的主要内容 2.1 随机信 ...查看


  • 初中物理电学经典难题整理

    • 初中物理--电学试题精选 1．(2015•盐亭县校级模拟)物理兴趣小组的同学独立组装成一盏台灯后,将其接入电路中．闭合开关时, 发现台灯不亮,同时其他用电器也停止工作,经检查发现保险丝熔断．更换新的保险丝后,将台灯单独接 入电路,又发现保险 ...查看


  • 初中物理电功率经典计算题50个

    • 初中物理电功率经典计算题100个(前50) 1.一根阻值为20Ω的电炉丝做成电热器,接上10V 电源能正常工作:用这个电热器与一个小灯泡串联接人该电源时,电热器消耗的功率是其额定功率的0.64倍.求: (1)电热器正常工作的电流和额定功率: ...查看


  • 发动机分为哪几类

    • 发动机分为哪几类 浏览次数:548次悬赏分:0 | 提问时间:2011-3-6 10:54 | 提问者:哭哉哀哉郭奉孝 目前应用于汽车的发动机主要有直列发动机,V 型发动机.W 型发动机.转子发动机几种类型.发动机根据点燃方式分为:压燃式和 ...查看


热门内容