FFT在功率谱密度计算中的应用
一、FFT算法理论依据和编程思想
FFT算法的基本思想：
考察DFT与IDFT的运算发现，利用以下两个特性可减少运算量：
Ⅰ)系数是一个周期函数，它的周期性和对称性可利用来改进运算，提高计算效率。如:
因此
利用这些周期性和对称性，DFT运算中有些项可合并；
Ⅱ）利用W N nk的周期性和对称性，把长度为N点的大点数的DFT运算分解为若干个小点数的DFT。因为DFT的计算量正比于N2，N小计算量也就小。FFT算是基于这样的基本思想发展起来的。它有多种形式，下面是按时间抽取的FFT（N点DFT运算的分解）先从一个特殊情况开始，假定N是2的整数次方，N=2M，M：正整数
1.将N点的DFT分解为两个N/2点的DFT：
首先将序列x（n）分解为两组，一组为偶数项，一组为奇数项
r=0,1，…，N/2-1
将DFT运算也相应分为两组：
其中X1（k）和X2（k）分别是x1（r）和x2（r）的N/2点DFT。
可见，一个N点的DFT可以分解为两个N/2点的DFT，这两个N/2点的DFT再按照上面（1）式合成为一个N点DFT，注意到，X1（k），X2（k）有N/2个点，即k=0，1，…，N/2-1，由（1）式得到X（k）只有N/2点，而实际上X（k）有N个点，即k=0，1，…，N-1，要用X1（k），X2（k）表示全部X（K）值，还必须应用系数w的周期性和对称性。
2.X(k)的(N/2)～N-1点表示：
由X(k)= X1（k）+w k N X2（k）, k=0,1,2,…,N/2-1
（2a）
得：
，
因为
,
且
同样
。
考虑到W N k对称性：。
故（2b）
（2a）式表示了X（k）前半部分k=0～N/2-1时的组成方式，（2b）式则表示了后半部分k=N/2～N-1时的组成方式。这两式所表示的运算过程可用一个称作蝶形的信号流图来表示。
3.蝶形信号流图：
如图1（a）所示，图中左面两支为输入，中间以一个小圆圈表示加、减运算，右上支为相加输出，右下支为相减输出，如果在某一支路上信号需要进行乘法运算，则在该支路上标以箭头，并将相乘的系数标在简头边，这样（2a），
（2b）所表示的运算，可用图1（b）所表示的“蝶形结”来表示。采用这种表示法，可将以上以讨论的分解过程用计算流图来表示。
图2.6所示为N=23=8的例子。通过这样分解以后，每一个N/2点DFT只需要
图2.6 N点DFT分解为2个N/2点DFT（N=8）
(N/2)2= N2/4次复数乘法运算，两个N/2点的DFT需要2(N/2)2= N2/2 次复乘，再加上将两个N/2点DFT合成为N点DFT时，在蝶形结前的N/2次复乘，共需要(N/2)2+N/2≈ N2/2次复乘，由此可见，经过这样的分解处理，运算量差不多节省了一倍。
4.将N/2点的DFT分解为两个N/4点的DFT：
既然这样分解是有效的，由于N=2M，N/2仍然是偶数，因此可对两个N/2点的DFT 再分别作进一步分解，例如对x1（r）和x2（r）可以再按其偶数部分及奇数部分分解为两个N/4点的DFT，既然这样分解是有效的，由于N=2M，N/2仍然是偶数，因此可对两个
N/2点的DFT再分别作进一步分解，例如对
x1（r）和x 2（r）可以再按其偶数部分及奇数
部分分解为两个N/4点的DFT，
l=0,1,…，N/4-1
而
同样X2(k)也可这样分解，并且将系数统一为，这样一个8点DFT就可分解为四个2点的DFT，如图2.7所示。
图2.7 N点DFT分解为4个N/4点的DFT(N=8)
5.2个点DFT的表示：
最后剩下的是2点DFT
，它可以用一个蝶形结表示，例如，x（0），x（4）所组成的2点DFT就可表示式：
这样，一个8点的完整的按时间抽取运算的流图如图2.8所示。
由于这样的方法每一步分解都是按输入时间序列是属于偶数还是奇数来抽取的，所以称为“按时间抽取法”或“时间抽取法”。
6.时间抽取法FFT运算特点：
（1）蝶形运算
对任何一2的整数幂N=2M，总是可以通过M次分解最后完全成为2点的DFT运算。这样的M次分解，就构成从x(n)到X(k)的M级运算过程。从上面的流图可看到，每一级
运算都由N/2个蝶形运算构成。因此每一级运算都需要
(N/2次复乘和N次复加(每个结
作加、减各一次)，这样，经过时间抽取后M级运算总共需要的运算：
复乘
复加N·M=Nlog2N
实际运算量与这个数字稍有出入，因为W 这几个系数实际上都不用乘法运算，因此在上面N=8的例子中，实际上只有两个系数W N1及W N3是需要乘法运算的。
用时间抽取法所需的计算量，不论是复乘还是复加都与Nlog2N成正比，而直接运算时则与N2成正比。例N=2048，N2=4194304，(N/2)log2N=11264，N2/[(N/2)log2N]=392.4倍。FFT显然要比直接法快得多。
(2）原位计算
当数据输入到存储器中以后，每一级运算的结果仍然储存在同一组存储器中，直到最后输出，中间无需其它存储器，这叫原位计算。
例如，N=8的FFT运算，输入x(0)，x(4)，x(2)，x(6)…，x(7)可分别存入A(1)，A(2)，…，A(8)这9个存储单元中，在第一级运算中，首先是存储单元A(1)，A(2)中x(0)，x(4)进入蝶形运算，x(0)，x(4)输入运算器后，其数值不再需要保存，因此蝶形运算的结果可仍然送回存储单元A(1)，A(2)中保存，然后A(3)，A(4)中x(2)，x(6)再进入蝶形运算，其结果再送回A(3)，A(4)，一直到算完A(7)，A(8)，则完成了第一级运算过程。第二级运算仍可采用这种原位的方式，但是进入蝶形结的组合关系不同，首先进入蝶形结的是A(1)、A(3)存储单元中的数据，运算结果仍可送回A(1)、A(3)保存，然后进入蝶形结的是A(2)、A(4)…，依此类推，每一级运算均可在原位进行，这种原位运算结构可节省存储单元，降低设备成本，还可节省找地址的时间。
（3）序数重排
对按时间抽取FFT的原位运算结构，当运算完毕时，这种结构存储单元A(1)、A(2)，…，A(8)中正好顺序存放着X(0)，X(1)，X(2)，…，X(7)，因此可直接按顺序输出，但这种原位运算的输入x(n)却不能按这种自然顺序存入存储单元中，而是按X(0)，X(4)，X(2)，
X(6)，…，X(7)的顺序存入存储单元，这种顺序看起来相当杂乱，然而它也是有规律的。当用二进制表示这个顺序时，它正好是“码位倒置”的顺序。例如，原来的自然顺序应是x(1)的地方，现在放着x(4)，用二进制码表示这一规律时，则是在x(0 0 1)处放着x(1 0 0)，x(0 1 1)处放着x(1 1 0)。即将自然顺序的二进制码位倒置过来，第一位码变成最末位码，这样倒置以后的顺序正是输入所需要的顺序。下表列出N=8时按码位倒置规律所得的顺序，其结果与按时间抽取FFT流图中的输入顺序是一致的。
表一码位倒置顺序
在实际运算中，一般直接将输入数据x（n）按码位倒置的顺序排好输入很不方便，总是先按自然顺序输入存储单元，然后再通过变址运算将自然顺序的存储转换成码位倒置顺
序的存储，然后进行FFT 的原位计算。目前有许多通用DSP 芯片支持这种码位倒置的寻址功能。
（4）蝶形类型随迭代次数成倍增加
观察8点FFT 的三次迭代运算
第一级迭代，只有一种类型的蝶形运算系数W 08
第二级迭代，有二种类型的蝶形运算系数W 08、W 28，参加运算的两个数据点间隔为2。 第三级迭代，有四类蝶形运算系数W 08、W 18、W 28、W 38，参加运算的两个数据点间隔为4。
所以，每次迭代的蝶形类型比上一次蝶代增加一倍，数据点间隔也增大一倍。 7.功率谱密度的计算
根据相关定理与维纳-辛钦关系式可得随机信号序列x(n)的功率谱密度
2
1()lim ()x N s k X k N
→∞= (1)
其估计值 2
1()()x S k X k N
∧
=
(2)
如果观察到序列x(n)的N 个值，即x(0), x(1),…x(N-1)，就可以通过FFT 直接求得X(k),再按式(2)求得S x (k),其计算过程如图2.9所示。
二、程序设计
本设计在主程序中分别调用了输入get_in()、倒序re_order()和蝶式运算butterfly(),功率谱密度计算power(),绘图on_draw 五个子函数。
1.程序框图 （1） 运算主程序框图
图2.9 周期图法计算功率谱密度流程
(2)、整序程序流程图
LH←N/2
j←LH
N1←N-2
i=1, N1
功率谱估计方法的比较 摘要： 本文归纳了信号处理中关键的一种分析方法, 即谱估计方法。概述了频谱估计中的周期图法、修正的协方差法和伯格递推法的原理，并且对此三种方法通过仿真做出了对比。 关键词：功率谱估计；AR 模型；参数 引言： 谱估计是指用已观测到的一定数量的样本数据估计一个平稳随机信号的谱。由于谱中包含了信号的很多频率信息，所以分析谱、对谱进行估计是信号处理的重要容。谱估计技术发展 渊源很长，它的应用领域十分广泛，遍及雷达、声纳、通信、地质勘探、天文、生物医学工程等众多领域，其容、方法都在不断更新，是一个具有强大生命力的研究领域。谱估计的理论和方法是伴随着随机信号统计量及其谱的发展而发展起来的，最早的谱估计方法是建 立在基于二阶统计量, 即自相关函数的功率谱估计的方法上。功率谱估计的方法经历了经典谱估计法和现代谱估计法两个研究历程，在过去及现在相当长一段时间里，功率谱估计一直占据着谱估计理论里的核心位置。经典谱估计也成为线性谱估计，包括BT 法、周期图法。现代谱估计法也称为非线性普估计，包括自相关法、修正的协方差法、伯格（Burg ）递推法、特征分解法等等。 原理： 经典谱估计方法计算简单，其主要特点是谱估计与任何模型参数无关，是一类非参数化的方法。它的主要问题是：由于假定信号的自相关函数在数据的观测区间以外等于零，因此估计出来的功率谱很难与信号的真实功率谱相匹配。在一般情况下，经典法的渐进性能无法给出实际功率谱的一个满意的近似，因而是一种低分辨率的谱估计方法。现代谱估计方法使用参数化的模型，他们统称为参数化功率谱估计，由于这类方法能够给出比经典法高得多的频率分辨率，故又称为高分辨率方法。下面分别介绍周期图法、修正的协方差法和伯格递推法。修正的协方差法和伯格递推法采用的模型均为AR 模型。 （1）周期图法 周期图法是先估计自相关函数, 然后进行傅里叶变换得到功率谱。假设随机信号x(n)只观测到一段样本数据，n=0, 1, 2, …, N-1。根据这一段样本数据估计自相关函数，如公式（1） 对（1）式进行傅里叶变换得到（2）式。 ∑--=+=1||0 *) ()(1 )(?m N n xx m n x n x N m r
功率谱定义 从确定性信号功率计算开始 ()()221 11lim lim 222T T T T T P x t dt X d T T ωωπ∞--∞→∞→∞==?? ()()21lim 2T T S X T ωω→∞= S(w)为功率谱密度，简称功率谱 则 ()12P S d ωωπ+∞-∞= ? 随机信号的功率谱密度 （1）样本功率谱与功率谱密度 ()()21,lim ,2X T T S X T ωξωξ→∞= 针对一个具体的样本而言，其是一个确定性的信号 (2) 随机信号的平均功率及平均功率谱密度 ()X X P E P ξ=???? 需要对具体的样本取概率均值才能计算出功率 ()()()21,lim ,2X X T T S E S E X T ωωξωξ→∞??==?????? 故功率谱密度是对所有概率取期望的反应。 （3）自相关函数与功率谱密度 ()()R S τω? （4）信号的自相关函数计算 分为确定信号和随机信号 确定信号 02002*0 1()lim ()()T T x T R x t x t dt T ττ-→∞=-? 周期信号 0202*0 1()()()T T x R x t x t dt T ττ-=-? 随机信号 *()[()()]x R E x t x t ττ=- 2 功率计算 （1）根据定义来计算
（2）周期信号如何计算 0cos()A t ω的计算 200()()1()[]2 A A s d T πσωωπσωωωω+∞-∞-++==?不好算因此放弃，但是应该可以类推得出结论 (3)自相关函数计算 0cos()A t ω的计算 /2 200/2 /222000/2201()cos()cos(())cos()cos(2)1[]2 cos()2 T T T T r A t t d T A A t d T A τωωτωωτωωτωωτ+-+-=-+-==?? 所以其功率谱为 200()2 A πσωωσωω（-）+（+） 0j t Ae ω的计算 0000/2()2/2 /22/2 21()1T j t j t T T j T j r A e e dt T A e dt T A e ωωτωτωτ τ+---+-===?? 总结：因此周期函数，首先转换成傅里叶级数，然后再通过自相关函数的定义计算自相关函数，得到其功率谱密度。
1 随机信号的经典谱估计方法 估计功率谱密度的平滑周期图是一种计算简单的经典方法。它的主要特点是与任 何模型参数无关，是一类非参数化方法[4]。它的主要问题是：由于假定信号的自相关函数在数据观测区以外等于零，因此估计出来的功率谱很难与信号的真实功率谱相匹配。在一般情况下，周期图的渐进性能无法给出实际功率谱的一个满意的近似，因而是一种低分辨率的谱估计方法。本章主要介绍了周期图法、相关法谱估计（BT ）、巴特利特(Bartlett)平均周期图的方法和Welch 法这四种方法。 2.1 周期图法 周期图法又称直接法。它是从随机信号x(n)中截取N 长的一段，把它视为能量有限x(n)真实功率谱)(jw x e S 的估计)(jw x e S 的抽样. 周期图这一概念早在1899年就提出了，但由于点数Ｎ一般比较大，该方法的计算量过大而在当时无法使用。只是1965年FFT 出现后，此法才变成谱估计的一个常用方法。周期图法[5]包含了下列两条假设： 1．认为随机序列是广义平稳且各态遍历的，可以用其一个样本x(n)中的一段 )(n x N 来估计该随机序列的功率谱。这当然必然带来误差。 2．由于对)(n x N 采用DFT ，就默认)(n x N 在时域是周期的，以及)(k x N 在频域是周期的。这种方法把随机序列样本x(n)看成是截得一段)(n x N 的周期延拓，这也就是周期图法这个名字的来历。与相关法相比，相关法在求相关函数)(m R x 时将 )(n x N 以外是数据全都看成零，因此相关法认为除)(n x N 外 x(n)是全零序列，这种处 理方法显然与周期图法不一样。 但是，当相关法被引入基于FFT 的快速相关后，相关法和周期图法开始融合。通过比较我们发现：如果相关法中M=N ，不加延迟窗，那么就和补充（N-1）个零的周期图法一样了。简单地可以这样说：周期图法是M=N 时相关法的特例。因此相关法和周期图法可结合使用。 2.2 相关法谱估计（BT ）法
t=0:0.0001:0.1; %时间间隔为0.0001，说明采样频率为10000Hz x=square(2*pi*1000*t); %产生基频为1000Hz的方波信号 n=randn(size(t)); %白噪声 f=x+n; %在信号中加入白噪声 figure(1); subplot(2,1,1); plot(f); %画出原始信号的波形图 ylabel('幅值(V)'); xlabel('时间(s)'); title('原始信号'); y=fft(f,1000); %对原始信号进行离散傅里叶变换，参加DFT采样点的个数为1000 subplot(2,1,2); m=abs(y); f1=(0:length(y)/2-1)'*10000/length(y);%计算变换后不同点对应的幅值plot(f1,m(1:length(y)/2)); ylabel('幅值的模'); xlabel('时间(s)'); title('原始信号傅里叶变换'); %用周期图法估计功率谱密度 p=y.*conj(y)/1000; %计算功率谱密度 ff=10000*(0:499)/1000; %计算变换后不同点对应的频率值 figure(2); plot(ff,p(1:500)); ylabel('幅值'); xlabel('频率(Hz)'); title('功率谱密度（周期图法）'); 功率谱估计在现代信号处理中是一个很重要的课题，涉及的问题很多。在这里，结合matlab，我做一个粗略介绍。功率谱估计可以分为经典谱估计方法与现代谱估计方法。经典谱估计中最简单的就是周期图法，又分为直接法与间接法。直接法先取N点数据的傅里叶变换（即频谱），然后取频谱与其共轭的乘积，就得到功率谱的估计；间接法先计
功率谱估计仿真实验 选题条件：对于给定的一个信号()()()t t f t f t x ?ππ++=212sin 2)2sin(，其中1f =50Hz ， 2f =100Hz ，()t ?为白噪声，采样频率Fs 为1000Hz ，对其进行功率谱估 计。 仿真目标：采用多种方法对该指定信号进行功率谱估计，计算其功率谱密度，比较 各种估计方法的优劣。 设计思路：本仿真实验采用经典谱估计中的周期图法对给定信号进行谱估计。但是 由于其自身的缺陷，使得频率分辨率较低。为了不断满足需要，找到恰 当的估计法，实验使依次使用了周期图法的改进型方法如分段周期图法、 窗函数法以及修正的周期图法进行功率谱估计，对四种方法得出的谱估 计波形进行比较分析，得出估计效果最好的基于周期图法的谱估计方法。 仿真指标：频率分辨率、估计量的方差、频谱光滑度 平台说明：本实验采用MATLAB7.0仿真软件，基于WINDOWS-XP 系统。Matlab 是 一个集数值分析、矩阵运算、信号处理和图形显示于一体的工程分析处理软件。它提供的部分算法函数为功率谱估计提供了一条可行的方便途径，如PSD 和CSD 可以自动实现Welch 法估计，而不需要自己编程。但是较为有限，大部分需要自己编写相应的M 文件来实现。 实现方法： 一、周期图法 周期图法是直接将信号的采样数据()n x 进行傅立叶变换求功率谱密度估计。假设有限长随机信号序列()n x ，将它的功率谱按定义写出如下： ()()??? ?????+=∑-=-∞→2121lim N N n n j N j xx e n x N E e P ωω 如果忽略上式中求统计平均的运算，观测数据为：()n x 10-≤≤N n ，便得到了周期图法的定义： ()()2 10 ^ 1n j N n j xx e n x N e P ωω--=∑=， 式中的绝对值符号内的部分可以用FFT 计算，这样就可得到周期图法的计算框图如下所示： () ω j xx e ^ 图1 周期图法计算功率谱框图 采用周期图法时，可以分取不同的信号长度256、512和1024，分别进行功率谱
功率谱估计及其MATLAB仿真 詹红艳 （201121070630控制理论与控制工程） 摘要:从介绍功率谱的估计原理入手分析了经典谱估计和现代谱估计两类估计方法的原理、各自特点及在Matlab中的实现方法。 关键词:功率谱估计;周期图法;AR参数法;Matlab Power Spectrum Density Estimation and the simulation in Matlab Zhan Hongyan （201121070630Control theory and control engineering） Abstract:Mainly introduces the principles of classical PSD estimation and modern PSD estimation,discusses the characteristics of the methods of realization in Matlab.Moreover,It gives an example of each part in realization using Matlab functions. Keywords:PSDPstimation,Periodogram method,AR Parameter method,Matlab 1引言 现代信号分析中,对于常见的具有各态历经的平稳随机信号,不可能用清楚的数学关系式来描述,但可以利用给定的N个样本数据估计一个平稳随机信号的功率谱密度叫做功率谱估计(PSD)。它是数字信号处理的重要研究内容之一。功率谱估计可以分为经典功率谱估计(非参数估计)和现代功率谱估计(参数估计)。 功率谱估计在实际工程中有重要应用价值,如在语音信号识别、雷达杂波分析、波达方向估计、地震勘探信号处理、水声信号处理、系统辨识中非线性系统识别、物理光学中透镜干涉、流体力学的内波分析、太阳黑子活动周期研究等许多领域,发挥了重要作用。 Matlab是MathWorks公司于1982年推出的一套高性能的数值计算和可视化软件,人称矩 阵实验室,它集数值分析、矩阵运算、信号处理和图形显示于一体,构成了一个方便的、界面友好的用户环境,成为目前极为流行的工程数学分析软件。也为数字信号处理进行理论学习、工程设计分析提供了相当便捷的途径。本文的仿真实验中,全部在Matlab6.5环境下调试通过;随机序列由频率不同的正弦信号加高斯白噪声组成。 2经典功率谱估计 经典功率谱估计是将数据工作区外的未知数据假设为零,相当于数据加窗。经典功率谱估计方法分为:相关函数法(BT法)、周期图法以及两种改进的周期图估计法即平均周期图法和平滑平均周期图法,其中周期图法应用较多,具有代表性。 1.1相关函数法(BT法) 该方法先由序列x(n)估计出自相关函数R(n),然后对R(n)进行傅立叶变换,便得到x(n)的功率谱估计。当延迟与数据长度相比很小时,可以有良好的估计精度。 Matlab代码示例1: Fs=500;%采样频率 n=0:1/Fs:1;
matlab实现功率谱密度分析psd及详细解说 功率谱密度幅值的具体含义？？ 求信号功率谱时候用下面的不同方法，功率谱密度的幅值大小相差很大！ 我的问题是，计算具体信号时，到底应该以什么准则决定该选用什么方法啊？ 功率谱密度的幅植的具体意义是什么？？下面是一些不同方法计算同一信号的matlab 程序！欢迎大家给点建议！ 直接法： 直接法又称周期图法，它是把随机序列x(n)的N个观测数据视为一能量有限的序列，直接计算x(n)的离散傅立叶变换，得X(k)，然后再取其幅值的平方，并除以N，作为序列x(n)真实功率谱的估计。 Matlab代码示例： clear; Fs=1000; %采样频率 n=0:1/Fs:1; %产生含有噪声的序列 xn=cos(2*pi*40*n)+3*cos(2*pi*100*n)+randn(size(n)); window=boxcar(length(xn)); %矩形窗 nfft=1024; [Pxx,f]=periodogram(xn,window,nfft,Fs); %直接法 plot(f,10*log10(Pxx)); 间接法： 间接法先由序列x(n)估计出自相关函数R(n)，然后对R(n)进行傅立叶变换，便得到x(n)的功率谱估计。 Matlab代码示例： clear; Fs=1000; %采样频率 n=0:1/Fs:1; %产生含有噪声的序列 xn=cos(2*pi*40*n)+3*cos(2*pi*100*n)+randn(size(n)); nfft=1024; cxn=xcorr(xn,'unbiased'); %计算序列的自相关函数 CXk=fft(cxn,nfft); Pxx=abs(CXk);
F F T在功率谱密度计算 中的应用 集团标准化工作小组 #Q8QGGQT-GX8G08Q8-GNQGJ8-MHHGN#
FFT在功率谱密度计算中的应用 一、FFT算法理论依据和编程思想 FFT算法的基本思想： 考察DFT与IDFT的运算发现，利用以下两个特性可减少运算量： Ⅰ)系数是一个周期函数，它的周期性和对称性可利用来改进运算，提高计算效率。如: 因此 利用这些周期性和对称性，DFT运算中有些项可合并； Ⅱ）利用W N nk的周期性和对称性，把长度为N点的大点数的DFT运算分解为若干个小点数的DFT。因为DFT的计算量正比于N2，N小计算量也就小。 FFT算法正是基于这样的基本思想发展起来的。它有多种形式，下面是按时间抽取的FFT（N点DFT运算的分解） 先从一个特殊情况开始，假定N是2的整数次方，N=2M，M：正整数 1.将N点的DFT分解为两个N/2点的DFT： 首先将序列x（n）分解为两组，一组为偶数项，一组为奇数项 r=0,1，…，N/2-1 将DFT运算也相应分为两组： 其中X 1（k）和X 2 （k）分别是x 1 （r）和x 2 （r）的N/2点DFT。 可见，一个N点的DFT可以分解为两个N/2点的DFT，这两个N/2点的DFT再按照上面 （1）式合成为一个N点DFT，注意到，X 1（k），X 2 （k）有N/2个点，即k=0，1，…， N/2-1，由（1）式得到X（k）只有N/2点，而实际上X（k）有N个点，即k=0， 1，…，N-1，要用X 1（k），X 2 （k）表示全部X（K）值，还必须应用系数w的周期性和 对称性。 (k)的(N/2)～N-1点表示： 由X(k)= X 1（k）+w k N X 2 （k）, k=0,1,2,…,N/2-1
1.基本方法 周期图法是直接将信号的采样数据x(n)进行Fourier变换求取功率谱密度估计的方法。假定有限长随机信号序列为x(n)。它的Fourier变换和功率谱密度估计存在下面的关系： 式中，N为随机信号序列x(n)的长度。在离散的频率点f=kΔf,有： 其中，FFT[x(n)]为对序列x(n)的Fourier变换，由于FFT[x(n)]的周期为N，求得的功率谱估计以N为周期，因此这种方法称为周期图法。下面用例子说明如何采用这种方法进行功率谱 用有限长样本序列的Fourier变换来表示随机序列的功率谱，只是一种估计或近似，不可避免存在误差。为了减少误差，使功率谱估计更加平滑，可采用分段平均周期图法（Bartlett法）、加窗平均周期图法（Welch 法）等方法加以改进。 2. 分段平均周期图法（Bartlett法） 将信号序列x(n)，n=0,1,…,N-1，分成互不重叠的P个小段，每小段由m个采样值，则P*m=N。对每个小段信号序列进行功率谱估计，然后再取平均作为整个序列x(n)的功率谱估计。 平均周期图法还可以对信号x(n)进行重叠分段，如按2:1重叠分段，即前一段信号和后一段信号有一半是重叠的。对每一小段信号序列进行功率谱估计，然后再取平均值作为整个序列x(n)的功率谱估计。这两种方法都称为平均周期图法，一般后者比前者好。程序运行结果为图9-5，上图采用不重叠分段法的功率谱估计，下图为2:1重叠分段的功率谱估计，可见后者估计曲线较为平滑。与上例比较，平均周期图法功率谱估计具有明显效果（涨落曲线靠近0dB）。 3.加窗平均周期图法 加窗平均周期图法是对分段平均周期图法的改进。在信号序列x(n)分段后，用非矩形窗口对每一小段信号序列进行预处理，再采用前述分段平均周期图法进行整个信号序列x(n)的功率谱估计。由窗函数的基本知识（第7章）可知，采用合适的非矩形窗口对信号进行处理可减小“频谱泄露”，同时可增加频峰的宽度，从而提高频谱分辨率。 其中上图采用无重叠数据分段的加窗平均周期图法进行功率谱估计，而下图采用重叠数据分段的加窗平均周期图法进行功率谱估计，显然后者是更佳的，信号谱峰加宽，而噪声谱均在0dB附近，更为平坦（注意采用无重叠数据分段噪声的最大的下降分贝数大于5dB，而重叠数据分段周期图法噪声的最大下降分贝数小于5dB）。 4. Welch法估计及其MATLAB函数 Welch功率谱密度就是用改进的平均周期图法来求取随机信号的功率谱密度估计的。Welch 法采用信号重叠分段、加窗函数和FFT算法等计算一个信号序列的自功率谱估计(PSD如上例中的下半部分的求法)和两个信号序列的互功率谱估计（CSD）。 MATLAB信号处理工具箱函数提供了专门的函数PSD和CSD自动实现Welch法估计，而不需要自己编程。 （1）函数psd利用Welch法估计一个信号自功率谱密度，函数调用格式为： [Pxx[,f]]=psd(x[,Nfft,Fs,window,Noverlap,’dflag’]) 式中，x为信号序列；Nfft为采用的FFT长度。这一值决定了功率谱估计速度，当Nfft采用2的幂时，程序采用快速算法；Fs为采样频率；Window定义窗函数和x分段序列的长度。窗函数长度必须小于或等于Nfft，否则会给出错误信息；Noverlap为分段序列重叠的采样
功率谱密度的估计 原始波=余弦波+白噪声 这个实验采用了两个输入，一个是白噪声，一个是有用信号和噪声信号作为输入时，他们的功率谱密度的仿真图像，并将他们进行对比。 平稳随机信号的功率谱密度（PSD ）是相关序列的离散傅里叶变换： ()()jw m XX x P w r m e ∞ --∞=∑ 采用间接法计算噪声信号的功率谱。 间接法，又称自相关法或者BT 法，在1985年由布莱克曼与图基首先开拓。间接法的理论基础是维纳-辛钦定理。他是由N 个观察值x(0),x(1),……,x(N-1),估计出自相关函数R （m ）,然后再求R （m ）的傅里叶变换作为功率谱密度的估计。 ()(),||1M jw jw m N m M S e R m e M N -=-=<=-∑ clear all; randn('state',0) NFFT=1024; %采样点数 Fs=1000; %取样频率（单位为Hz ） t=0:1/Fs:.2;
y1=cos(t*20*pi); %余弦序列 figure(1) plot(t,y1); ylabel('余弦序列'); grid on; %余弦序列的图像： %白噪声 m=(0:NFFT-1)/Fs; y=0.1*randn(size(m)); %产生高斯白噪声。 figure(2); plot(m,y); title('白噪声波形'); grid on;
%白噪声的自相关函数 [cory,lags]=xcorr(y,200,'unbiased'); %计算白噪声的自相关函数 figure(3) plot(lags,cory); %自相关函数(无偏差的)，其中，cory为要求的自相关函数，lag为自相关函数的长度。 title('白噪声相关函数'); grid on;
清华大学 《数字信号处理》期末作业 2013 年 1 月
第一题掌握去噪的方法 1.1 题目描述 MATLAB 中的数据文件noisdopp 含有噪声，该数据的抽样频率未知。调出该数据，用你学过的滤波方法和奇异值分解的方法对其去噪。要求：1．尽可能多地去除噪声，而又不损害原信号； 2．给出你去噪的原理与方法；给出说明去噪效果的方法或指标； 3．形成报告时应包含上述内容及必要的图形，并附上原程序。 1.2 信号特性分析 MATLAB所给noisdopp信号极其频域特征如图1.1、图1.2。 图1.1含有噪声的noisdopp信号
图1.2 noisdopp 信号频域特性 其中横坐标f 采用归一化频率，即未知抽样频率Fs 对应2（与滤波器设计时参数一致）。信号基本特性是一个幅值和频率逐渐增加的正弦信号叠加噪声，噪声为均匀的近似白噪声，没有周期等特点。 因为噪声信号能量在全频带均匀分布，滤波器截止频率过低则信号损失大，过高则噪声抑制小，认为频谱中含有毛刺较多的部分即为信噪比较小的部分，滤除这部分可以达到较好的滤波效果。 先给定去噪效果的评定指标。信号开始阶段频率较高（如图1.3，红圈为信号值），一周期内采样点4~5个，即信号归一化频率达到0.4~0.5（Fs=2），难以从频域将这部分信号同噪声分离，滤波后信号损失较大，故对前128点用信噪比考察其滤波效果，定义： 2 2 () 10lg (()())k k x k SNR y k x k =-∑∑ 其中，()x k 为原nosidopp 信号，()y k 为滤波后信号。SNR 越大表示滤除部分能力越小，可以反映滤波后信号对原信号的跟踪能力，对前128点主要考察SNR ，越大滤波器性能越好。
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振幅谱平方的平均值，并转化为dB f=(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振幅谱平方的平均值，并转化为dB f=(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振幅谱平方的平均值，并转化为dB f=(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;
功率谱估计性能分析及Matlab 仿真 1 引言 随机信号在时域上是无限长的，在测量样本上也是无穷多的，因此随机信号的能量是无限的，应该用功率信号来描述。然而，功率信号不满足傅里叶变换的狄里克雷绝对可积的条件，因此严格意义上随机信号的傅里叶变换是不存在的。因此，要实现随机信号的频域分析，不能简单从频谱的概念出发进行研究，而是功率谱[1]。 信号的功率谱密度描述随机信号的功率在频域随频率的分布。利用给定的 N 个样本数据估计一个平稳随机信号的功率谱密度叫做谱估计。谱估计方法分为两大类：经典谱估计和现代谱估计。经典功率谱估计如周期图法、自相关法等，其主要缺陷是描述功率谱波动的数字特征方差性能较差，频率分辨率低。方差性能差的原因是无法获得按功率谱密度定义中求均值和求极限的运算[2]。分辨率低的原因是在周期图法中，假定延迟窗以外的自相关函数全为0。这是不符合实际情况的，因而产生了较差的频率分辨率。而现代谱估计的目标都是旨在改善谱估计的分辨率，如自相关法和Burg 法等。 2 经典功率谱估计 经典功率谱估计是截取较长的数据链中的一段作为工作区，而工作区之外的数据假设为0，这样就相当将数据加一窗函数，根据截取的N 个样本数据估计出其功率谱[1]。 2.1 周期图法( Periodogram ) Schuster 首先提出周期图法。周期图法是根据各态历经的随机过程功率谱的定义进行的谱估计。 取平稳随机信号()x n 的有限个观察值(0),(1),...,(1)x x x n -，求出其傅里叶变换 1 ()()N j j n N n X e x n e ω ω---==∑ 然后进行谱估计
Classical Power Spectrum Estimation Abstract With the increasing need of spectrum, various computational methods and algorithms have been proposed in the literature. Keeping these views and facts of spectrum shaping capability by FRFT based windows we have proposed a closed form solution for Bartlett window in fractional domain. This may be useful for analysis of different upcoming generations of mobile communication in a better way which are based on OFDM technique. Moreover, it is useful for real-time processing of non-stationary signals. As per our best knowledge the closed form solution mentioned in this paper have not been reported in the literature till date.This paper focuses on classical period spectral estimation and moderu spectral estimation based on Burg algorithm. By comparing various algorithms in computational complexity and resolution, Burg algorithm was used to signal processing finally. Experimental and simulation results indicated that digital signal processing system would meet system requirements for measurement accuracy. Keywords periodogram spectral estimation ; Burg algorithm I. INTRODUCTION When we expand the frequency response of any digital filter by means of Fourier series, we get impulse response of the digital filter in the form of coefficients of the Fourier series. But the resultant filter is unrealizable and also its impulse response in infinite in duration. If we directly truncate this series to a finite number of points we have to face with well known Gibbs phenomenon, so we modify the Fourier coefficients by
FFT在功率谱密度计算中的应用 一、FFT算法理论依据和编程思想 FFT算法的基本思想： 考察DFT与IDFT的运算发现，利用以下两个特性可减少运算量： Ⅰ)系数是一个周期函数，它的周期性和对称性可利用来改进运算，提高计算效率。如: 因此 利用这些周期性和对称性，DFT运算中有些项可合并； Ⅱ）利用W N nk的周期性和对称性，把长度为N点的大点数的DFT运算分解为若干个小点数的DFT。因为DFT的计算量正比于N2，N小计算量也就小。FFT算是基于这样的基本思想发展起来的。它有多种形式，下面是按时间抽取的FFT（N点DFT运算的分解）先从一个特殊情况开始，假定N是2的整数次方，N=2M，M：正整数 1.将N点的DFT分解为两个N/2点的DFT： 首先将序列x（n）分解为两组，一组为偶数项，一组为奇数项 r=0,1，…，N/2-1 将DFT运算也相应分为两组： 其中X1（k）和X2（k）分别是x1（r）和x2（r）的N/2点DFT。 可见，一个N点的DFT可以分解为两个N/2点的DFT，这两个N/2点的DFT再按照上面（1）式合成为一个N点DFT，注意到，X1（k），X2（k）有N/2个点，即k=0，1，…，N/2-1，由（1）式得到X（k）只有N/2点，而实际上X（k）有N个点，即k=0，1，…，N-1，要用X1（k），X2（k）表示全部X（K）值，还必须应用系数w的周期性和对称性。
2.X(k)的(N/2)～N-1点表示： 由X(k)= X1（k）+w k N X2（k）, k=0,1,2,…,N/2-1 （2a） 得： ， 因为 , 且 同样 。 考虑到W N k对称性：。 故（2b） （2a）式表示了X（k）前半部分k=0～N/2-1时的组成方式，（2b）式则表示了后半部分k=N/2～N-1时的组成方式。这两式所表示的运算过程可用一个称作蝶形的信号流图来表示。 3.蝶形信号流图：
功率谱密度的三种matlab 实现方法 一：实验目的： (1)掌握三种算法的概念、应用及特点； (2)了解谱估计在信号分析中的作用； (3) 能够利用burg 法对信号作谱估计，对信号的特点加以分析。 二;实验内容： （1）简单说明三种方法的原理。 （2）用三种方法编写程序，在matlab 中实现。 （3）将计算结果表示成图形的形式，给出三种情况的功率谱图。 （4）比较三种方法的特性。 （5）写出自己的心得体会。 三：实验原理： 1.周期图法： 周期图法又称直接法。它是从随机信号x(n)中截取N 长的一段，把它视为能量有限x(n)真实功率谱)(jw x e S 的估计)(jw x e S 的抽样. 认为随机序列是广义平稳且各态遍历的，可以用其一个样本x(n)中的一段)(n x N 来估计该随机序列的功率谱。这当然必然带来误差。由于对)(n x N 采用DFT ，就默认)(n x N 在时域是周期的，以及)(k x N 在频域是周期的。这种方法把随机序列样本x(n)看成是截得一段)(n x N 的周期延拓，这也就是周期图法这个名字的来历。
2.相关法（间接法): 这种方法以相关函数为媒介来计算功率谱，所以又叫间接法。这种方法的具体步骤是： 第一步：从无限长随机序列x(n)中截取长度N 的有限长序列列 )(n x N 第二步：由N 长序列)(n x N 求（2M-1）点的自相关函数)(m R x ∧ 序列。 )()(1)(1 m n x n x N m R N n N N x += ∑-=∧ (2-1) 这里，m=-(M-1)…,-1,0,1…,M-1，M N ，)(m R x 是双边序列，但是由自相关函数的偶对称性式，只要求出m=0,。。。,M-1的傅里叶变换，另一半也就知道了。 第三步：由相关函数的傅式变换求功率谱。即 jwm M M m X jw x e m R e S ----=∧∧ ∑= )()(1) 1( 以上过程中经历了两次截断，一次是将x(n)截成N 长，称为加数据窗，一次是将x(n)截成（2M-1）长，称为加延迟窗。因此所得的功率谱仅是近似值，也叫谱估计，式中的)(jw x e S 代表估值。一般取M< 功率谱估计浅谈 摘要：介绍了几种常用的经典功率谱估计与现代功率谱估计的方法原理，并利用Matlab对随机信号进行功率谱估计，对两种方法做出比较，分别给出其优缺点。关键词：功率谱；功率谱估计；经典功率谱估计；现代功率谱估计 前言 功率谱估计是从频率分析随机信号的一种方法,一般分成两大类:一类是经典谱估计;另一类是现代谱估计。由于经典谱估计中将数据工作区以外的未知数据假设为零，这相当于数据加窗，导致分辨率降低和谱估计不稳定。现代谱估计则不再简单地将观察区外的未知数据假设为零，而是先将信号的观测数据估计模型参数，按照求模型输出功率的方法估计信号功率谱，回避了数据观测区以外的数据假设问题。 周期图、自相关法及其改进方法(Welch)为经典(非参数)谱估计方法, 其以相关和傅里叶变换为基础,对于长数据记录较适用,但无法根本解决频率分辨率低和谱估计稳定性的问题,特别是在数据记录很短的情况下,这一问题尤其突出。以随机过程的参数模型为基础的现代参数法功率谱估计具有更高的频率分辨率和更好的适应性,可实现信号检测或信噪分离,对语音、声纳雷达、电磁波及地震波等信号处理具有重要意义,并广泛应用于通信、自动控制、地球物理等领域。在现代参数法功率谱估计方法中,比较有效且实用的是AR模型法，Burg谱估计法,现代谱估计避免了计算相关,对短数据具有更强的适应性,从而弥补了经典谱估计法的不足,但其也有一些自身的缺陷。 下面就给出这两类谱估计的简单原理介绍与方法实现。 经典谱估计法 经典法是基于传统的傅里叶变换。本文主要介绍一种方法：周期图法。 周期图法 由于对信号做功率谱估计，需要用计算机实现，如果是连续信号，则需要变换为离散信号。下面讨论离散随机信号序列的功率谱问题。 连续时间随机信号的功率谱密度与自相关函数是一对傅里叶变换对，即： 现代数字信号处理功率谱估计实验 12040007陈骁 一、实验条件 使用MATLAB 中数字信号处理工具的功率谱估计函数对信号进行估计。periodogram,pwelch,pyulear,pburg,pcov,pmcov,pmusic 函数分别对应课堂上讲的周期图法、welch 法、自相关法、协方差法、修正协方差法和burg 法。 本次实验用一个白噪声通过五阶的AR 过程产生一个有色噪声，功率为1.21。信号设计为：()()()1000sin 0.11000sin 0.110.2sin 0.4x t t t πππ=++。即在0.1π和0.11π的位置上有两个频率接近的大功率的信号，信噪比为：53.9dB ，在0.4π的位置上有一个小信号，信噪比为：-20dB 。产生的信号共计600个点。 用各种方法、每种方法用几种不同阶数进行估计。绘制估计得到的功率谱。分析各种方法的效果。 实际信号如图： 实际的功率谱如图： 二、实验结果 1、用周期图法得到的功率谱 2、welch方法 1、将数据分成8段 2、将数据分成3段 3、自相关法 1、20阶 2、50阶 4、burg法 1、20阶 2、50阶 5、music法 1、10阶 1、20阶 3、100阶 三、结果分析 1、普通的周期图法得到的功率谱确实方差性能很差，起伏十分剧烈。其可以将两个相邻的强信号分辨出来，分辨率不错。但是，频率在0.4处小信号，由于功率谱的起伏太剧烈，所以信号几乎无法被识别出来。 2、用Welch方法求出周期图，其段间重叠为默认值50%。分别将数据分为8段和3段。在将数据分成8段时，在频率为0.1附近的两个信号完全无法区分出来，在频率为0.4处的信号也几乎无法被检测到。而将数据分为3段时，两个间距较近的信号已经出现了两个峰值，说明其分辨率较8段时好，但比用普通的周期图法要差不少，在频率为0.4处，信号也有所表现。在功率谱的方差性能上，用8段的比用3段的强，又比不分段的强。 3、用AR模型的自相关法，分别取20阶、50阶、100阶得到周期图。可以看出来：阶数越高时功率谱的分辨率就会越高，在20阶时相邻的信号完全无法分辨，50阶时已经开始分离，100阶时已经分离的比较好了。而对于弱信号，用这个方 功率谱密度幅值的具体含义？？ 求信号功率谱时候用下面的不同方法，功率谱密度的幅值大小相差很大！ 我的问题是，计算具体信号时，到底应该以什么准则决定该选用什么方法啊？ 功率谱密度的幅植的具体意义是什么？？下面是一些不同方法计算同一信号的matlab 程序！欢迎大家给点建议！ 一、直接法： 直接法又称周期图法，它是把随机序列x(n)的N个观测数据视为一能量有限的序列，直接计算x(n)的离散傅立叶变换，得X(k)，然后再取其幅值的平方，并除以N，作为序列x(n)真实功率谱的估计。 Matlab代码示例： clear; Fs=1000; %采样频率 n=0:1/Fs:1; %产生含有噪声的序列 xn=cos(2*pi*40*n)+3*cos(2*pi*100*n)+randn(size(n)); window=boxcar(length(xn)); %矩形窗 nfft=1024; [Pxx,f]=periodogram(xn,window,nfft,Fs); %直接法 plot(f,10*log10(Pxx)); 二、间接法： 间接法先由序列x(n)估计出自相关函数R(n)，然后对R(n)进行傅立叶变换，便得到x(n)的功率谱估计。 Matlab代码示例： clear; Fs=1000; %采样频率 n=0:1/Fs:1; %产生含有噪声的序列 xn=cos(2*pi*40*n)+3*cos(2*pi*100*n)+randn(size(n)); nfft=1024; cxn=xcorr(xn,'unbiased'); %计算序列的自相关函数 CXk=fft(cxn,nfft); Pxx=abs(CXk); index=0:round(nfft/2-1); k=index*Fs/nfft; plot_Pxx=10*log10(Pxx(index+1)); plot(k,plot_Pxx); 三、改进的直接法： 对于直接法的功率谱估计，当数据长度N太大时，谱曲线起伏加剧，若N太小，谱的分辨率又不好，因此需要改进。功率谱估计浅谈汇总
数字信号处理——功率谱估计
计算功率谱密度
相关主题
文本预览