[转载]功率谱

原文地址：功率谱作者：wangjiufa1987

求信号功率谱时候用下面的不同方法，功率谱密度的幅值大小相差很大！
我的问题是，计算具体信号时，到底应该以什么准则决定该选用什么方法啊？
功率谱密度的幅植的具体意义是什么？？下面是一些不同方法计算同一信号的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太小，谱的分辨率又不好，因此需要改进。

1. Bartlett法

Bartlett平均周期图的方法是将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;
window=boxcar(length(n)); %矩形窗
noverlap=0; %数据无重叠
p=0.9; %置信概率

[Pxx,Pxxc]=psd(xn,nfft,Fs,window,noverlap,p);

index=0:round(nfft/2-1);
k=index*Fs/nfft;
plot_Pxx=10*log10(Pxx(index+1));
plot_Pxxc=10*log10(Pxxc(index+1));
figure(1)
plot(k,plot_Pxx);

pause;

figure(2)
plot(k,[plot_Pxx plot_Pxx-plot_Pxxc plot_Pxx+plot_Pxxc]);

2. Welch法

Welch法对Bartlett法进行了两方面的修正，一是选择适当的窗函数w(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;
window=boxcar(100); %矩形窗
window1=hamming(100); %海明窗
window2=blackman(100); %blackman窗
noverlap=20; %数据无重叠
range='half'; %频率间隔为[0 Fs/2]，只计算一半的频率

[Pxx,f]=pwelch(xn,window,noverlap,nfft,Fs,range);
[Pxx1,f]=pwelch(xn,window1,noverlap,nfft,Fs,range);
[Pxx2,f]=pwelch(xn,window2,noverlap,nfft,Fs,range);

plot_Pxx=10*log10(Pxx);
plot_Pxx1=10*log10(Pxx1);
plot_Pxx2=10*log10(Pxx2);

figure(1)
plot(f,plot_Pxx);

pause;

figure(2)
plot(f,plot_Pxx1);

pause;

figure(3)
plot(f,plot_Pxx2);

 

功率谱的数据都是相对值，他无法给出信号的实际绝对幅值，一般只要看峰值之间的比值正确就行了，当然这个问题可以通过做正规化处理解决

 

谢谢回答!不过具体原因,我还是不很明白啊?你能不能从原理上讲讲看啊!而且有时候所求信号的幅值意义是很大的!比如,地震信号的功率谱值,其幅值有一定的范围,而我求出来的值总是和文献的对不上,不知道具体选择求解方法时怎么处理啊??

 

xn=cos(2*pi*40*n)+3*cos(2*pi*100*n)
像这个信号的话,考虑单边谱的话,40Hz的幅值为1,100Hz的幅值为3,对应的功率谱值分别为1和9.
在用FFT做谱估计值时,应把FFT的结果取模后除以FFT的点数再乘以2,得到单边谱幅值,再平方后
就得到单边功率谱值

 

大家继续讨论啊!在地震信号处理中，常常统计功率谱幅值的变化规律，按大家的说法，就毫无意义了，因为采用不同方法，幅值差别很大，到底是怎么一回事啊？？

 

实际上，功率谱也可以求出真实幅值的，只要把求得的结果对
采样频率归一化，即按楼主的方法得到的相对值再乘以采样频率
即为真实结果！

 

用各种方法得到的功率谱，都是相对值，往往用分贝来表示。则是否能知道它的绝对值，这是可能的，但需要进行校正，要设定0分贝对应的数值。
我们测量的无非是加速度、速度和位移的功率谱，则它们0分贝的参考值分别是：
a0=1um/s^2
v0=1nm/s
d0=1pm
（摘自马大猷《声学手册》）。为了能求出绝对的功率谱密度数值，要从测量源头开始便要进行校正，也就是从传感器、放大器、直至进AD变换，每一环节都要校正，这样才能拿到正确的结果。

 

功率谱密度，单位为：unit^2/Hz代表单位频率上信号的能量，所以是密度谱，幅值代表频段内的有效值平方，计算时的步骤为
1 对每一分段数据进行FFT变换，并求的幅值谱
2 对幅值谱进行平方
3 将双边谱转化为单边谱
4 除以频率分辨率
举个例子：
幅值为1，频率为16Hz的正弦信号，使用1024Hz采样，2048点进行功率谱密度计算，频率分辨率为1024/2048=0.5Hz，求出的功率谱单边谱在第32根谱线处的值为1，解释为：信号FFT变换后得到的双边谱，幅值分别为0.5，平方后为0.25，转化为单边乘2为0.5，在除以频率分辨率为1。将1乘以0.5，正好为该信号有效值0.707的平方。

 

能不能解释一下为什么要转化为单边谱，而且最后除0.5是什么意思？谢谢

 

因为实数信号的双边幅值谱都是对称的，因此用单边谱就够了，这时候把负频率成分附加到相应的正频率成分，也就是在双边谱的基础上乘以２．

 

双边谱中包含负频率，在物理系统中是没有的；0.5为频率分辨率。

      引用:

原帖由 yangzp 于 2006-9-23 15:31 发表
功率谱密度，单位为：unit^2/Hz代表单位频率上信号的能量，所以是密度谱，幅值代表频段内的有效值平方，计算时的步骤为
1 对每一分段数据进行FFT变换，并求的幅值谱
2 对幅值谱进行平方
3 将双边谱转化为单边谱 ...

我周期图法得到的幅值为1，频率为16Hz的正弦信号的功率谱单边谱在第32根谱线处的值不为1呀？
能不能解释一下，用周期图法得到在某一个频率下的功率谱与信号的幅值有什么关系？谢谢！！

 

以上针对楼上的代码

问题1:
[Pxx,f]=pwelch(xn,window,noverlap,nfft,Fs,range);
[plot_Pxx=10*log10(Pxx);
我一直想知道功率谱计算完毕.为什么要取对数?
问题2:
[Pxx,Pxxc]=psd(xn,nfft,Fs,window,noverlap,p);
plot_Pxxc=10*log10(Pxxc(index+1));
问什么要将index索引+1呢?
这样+1后,谱线不就变成511条而不是512条了吗?
而用[Pxx,f]=pwelch(xn,window,noverlap,nfft,Fs,range); 方法没有对数组索引进行+1移位呢.

 

问题1：Pxx是功率谱，它的表示方法可以用线性，也可以用对数。因为在功率谱中变化的范围可能跨越几个数量级，用线性标度表示时，看不出这个特性，而用对数标度表示便能更清楚地表示出来。
问题2，index的定义为：
index=0:round(nfft/2-1);
共有512个数，在数组Pxxc中下标不能从0开始，必须从1开始，这便是为什么要设成Pxxc(index+1)，一样有512个数。

 

上面几位讲的非常精彩，但还有一个问题，就是：
      原始数据是d()，那么a=fft(d)得到的是它的傅立叶变换，在求功率谱的时候，应该是先求a的模，然后除以fft的点数，再下面问题来了，前面几位产生了不同的说法：
      1、yangzj说得是先乘以2，得到单边的，然后再平方，再除以频率分辨率，得到的就是单边功率谱。
      他的例子是：xn=cos(2*pi*40*n)+3*cos(2*pi*100*n)像这个信号的话,考虑单边谱的话,40Hz的幅值为1,100Hz的幅值为3,对应的功率谱值分别为1和9.在用FFT做谱估计值时,应把FFT的结果取模后除以FFT的点数再乘以2,得到单边谱幅值,再平方后就得到单边功率谱值。
      2、yangzp说得是先平方，然后乘以2，再除以频率分辨率，得到的就是单边功率谱。
      他的例子是：幅值为1，频率为16Hz的正弦信号，使用1024Hz采样，2048点进行功率谱密度计算，频率分辨率为1024/2048=0.5Hz，求出的功率谱单边谱在第32根谱线处的值为1，解释为：信号FFT变换后得到的双边谱，幅值分别为0.5，平方后为0.25，转化为单边乘2为0.5，在除以频率分辨率为1。将1乘以0.5，正好为该信号有效值0.707的平方。
       大家看，他们两个说的明显的差一个二倍，到底那个是对的呢？

 

1)求信号功率谱时候用下面的不同方法，功率谱密度的幅值大小相差很大！
  matlab提供的功率谱求法肯定是没有错的，关键在于各个参数的选取，以及对每种求解方法的了解！
平均周期图法和其他方法求出的结果，参数条件取得一样的话，可以得到完全相同的结果。
2）我的问题是，计算具体信号时，到底应该以什么准则决定该选用什么方法啊？
   直接采用平均周期图法求功率谱时，功率普形状呈锯齿形，谱峰点的准确位置不大好定。于是可以采用其他的方法对谱进行平滑操作，平滑化，仅仅是为了使图形光滑，并不会使得波的本质受到歪曲和畸变。反过来说，由于不纯的东西去掉了，本质的东西必然会更加显示出来！平滑化的程度可以根据所分析的信号，选择合理的窗函数和带宽！可以采用逐步观察的方法进行考察！
3）功率谱密度的幅植的具体意义是什么？？
   对于地震动信号（我研究的范畴），对于功率谱的单位，如数据为加速度时，为米（平方）/秒（三次方）；若数据为速度或者位移时，他的单位为米（平方）/秒，米（平方）*秒；他和物理学中的功率（单位时间内所作的功，单位为瓦特）概念是不一样的。因此在地震波的情况下，所谓功率谱，不能用功率这一物理量来理解。
    对于一个问题的理解，关键靠自己多实践，找相关的资料学习，这样可以加深理解

 

我现在通过试验得到离散化的峰值谱曲线，按上面讨论的提示，是不是只要先进行平方，然后除以频率分辨率就可以了？我们试验测试的都是单边谱，是否还需要×2。谢谢大家解答。

 

你肯定理解错了单边谱和双边谱的概念了。单边谱就是双边谱乘以2，所以你得到的如果是单边谱的话，那么应该在平方然后除以频率分辨率后，再除以2

 

     引用:

原帖由 yangzp 于 2006-9-23 15:31 发表
功率谱密度，单位为：unit^2/Hz代表单位频率上信号的能量，所以是密度谱，幅值代表频段内的有效值平方，计算时的步骤为
1 对每一分段数据进行FFT变换，并求的幅值谱
2 对幅值谱进行平方
3 将双边谱转化为单边谱
...

按[Pxx,f]=pwelch(xn,window,noverlap,nfft,fs,range);求得的Pxx是不是等于1/N*|X(k)|.^2?即求到了上面步骤2
如果我需要的是单边谱，则为Pxx*2，即为上面步骤3
这样计算之后，是否还需再除以采样频率？

 

其实matlab的结果是已经除过采样频率的，见pwelch的帮助文件，matlab求的是S(exp(jw))/Fs，应该说结果该不该乘以Fs

 

自相关函数DFT计算

这里关键是离散时, 间接法如何计算自相关函数和功率谱是一对Fourier变换.

如楼主程序中,取xn=1024点,直接FFT得功率谱是1024点(为直观,简单,取Fs=nfft=1024;)

xn自相关后是2047点,但楼主程序中对其作1024点FFT,得功率谱1024点. (即取Xn的前1024点), 问题是自相关后是2047点,什么得1024谱线?

一般数据DFT中的n=0:N-1;  在计算自相关函数DFT时, n=-N+1:N-1, 这样2047个样点得谱线是1024个.

下面的程序用直接法和间接法结果相同了, 为简单,取Fs=N=1024;直接法fft后计功率谱,

close all;clc;clear all;
Fs=1024; %采样频率
n=0:1/Fs:1-1/Fs;
%产生含有噪声的序列
xn=cos(2*pi*40*n)+3*cos(2*pi*100*n)+randn(size(n));
nfft=1024;
Pxx=fft(xn).*conj(fft(xn));%直接法
plot(10*log10(Pxx),'b');
hold
Fs=1024; %采样频率
n=0:1/Fs:1-1/Fs;
%产生含有噪声的序列
xn=cos(2*pi*40*n)+3*cos(2*pi*100*n)+randn(size(n));
nfft=1024;
cxn=xcorr(xn);%计算序列的自相关函数
cxn=cxn(nfft:end)+[0 cxn(1:nfft-1)];
CXk=fft(cxn,nfft);
Pxx=abs(CXk);
plot_Pxx=10*log10(Pxx);
plot(plot_Pxx,'r');

图中直接法功率谱(兰)  间接法功率谱(红)

 

经典谱估计方法的MATLAB分析
姚武川　姚天任

 

功率谱密度谱是一种概率统计方法，是对随机变量均方值的量度。一般用于随机振动分析，连续瞬态响应只能通过概率分布函数进行描述，即出现某水平响应所对应的概率。

功率谱密度是结构在随机动态载荷激励下响应的统计结果，是一条功率谱密度值—频率值的关系曲线，其中功率谱密度可以是位移功率谱密度、速度功率谱密度、加速度功率谱密度、力功率谱密度等形式。数学上，功率谱密度值—频率值的关系曲线下的面积就是方差，即响应标准偏差的平方值。

谱是个很不严格 的东西，常常指信号的Fourier变换，是一个时间平均（time average）概念功率谱的概念是针对功率有限信号的(能量有限信号可用能量谱分析)，所表现的是单位频带内信号功率随频率的变换情况。保留频谱的幅度信息，但是丢掉了相位信息，所以频谱不同的信号其功率谱是可能相同的。有两个重要区别： 1。功率谱是随机过程的统计平均概念，平稳随机过程的功率谱是一个确定函数；而频谱是随机过程样本的Fourier变换，对于一个随机过程而言，频谱也是一个“随机过程”。（随机的频域序列） 2。功率概念和幅度概念的差别。此外，只能对宽平稳的各态历经的二阶矩过程谈功率谱，其存在性取决于二阶局是否存在并且二阶矩的Fourier变换收敛；而频谱的存在性仅仅取决于该随机过程的该样本的Fourier变换是否收敛。

频谱分析（也称频率分析），是对动态信号在频率域内进行分析，分析的结果是以频率为坐标的各种物理量的谱线和曲线，可得到各种幅值以频率为变量的频谱函数F(ω)。频谱分析中可求得幅值谱、相位谱、功率谱和各种谱密度等等。频谱分析过程较为复杂，它是以傅里叶级数和傅里叶积分为基础的。

功率谱是个什么概念？它有单位吗?

随机信号是时域无限信号，不具备可积分条件，因此不能直接进行傅氏变换。一般用具有统计特性的功率谱来作为谱分析的依据。功率谱与自相关函数是一个傅氏变换对。功率谱具有单位频率的平均功率量纲。所以标准叫法是功率谱密度。通过功率谱密度函数，可以看出随机信号的能量随着频率的分布情况。像白噪声就是平行于 w轴，在w轴上方的一条直线。

功率谱密度，从名字分解来看就是说，观察对象是功率，观察域是谱域，通常指频域，密度，就是指观察对象在观察域上的分布情况。一般我们讲的功率谱密度都是针对平稳随机过程的，由于平稳随机过程的样本函数一般不是绝对可积的，因此不能直接对它进行傅立叶分析。可以有三种办法来重新定义谱密度，来克服上述困难。一是用相关函数的傅立叶变换来定义谱密度；二是用随机过程的有限时间傅立叶变换来定义谱密度；三是用平稳随机过程的谱分解来定义谱密度。三种定义方式对应于不同的用处，首先第一种方式前提是平稳随机过程不包含周期分量并且均值为零，这样才能保证相关函数在时差趋向于无穷时衰减，所以lonelystar说的不全对，光靠相关函数解决不了许多问题，要求太严格了；对于第二种方式，虽然一个平稳随机过程在无限时间上不能进行傅立叶变换，但是对于有限区间，傅立叶变换总是存在的，可以先架构有限时间区间上的变换，在对时间区间取极限，这个定义方式就是当前快速傅立叶变换（FFT）估计谱密度的依据；第三种方式是根据维纳的广义谐和分析理论：Generalized harmonic analysis, Acta Math, 55(1930),117-258,利用傅立叶-斯蒂吉斯积分，对均方连续的零均值平稳随机过程进行重构，在依靠正交性来建立的。

另外，对于非平稳随机过程，也有三种谱密度建立方法，由于字数限制，功率谱密度的单位是G的平方/频率。就是就是函数幅值的均方根值与频率之比。是对随机振动进行分析的重要参数。

功率谱密度的国际单位是什么?

如果是加速度功率谱密度，加速度的单位是m/s^2,
那么，加速度功率谱密度的单位就是(m/s^2)^2/Hz,
而Hz的单位是1/s,经过换算得到加速度功率谱密度的单位是m^2/s^3.

同理，如果是位移功率谱密度，它的单位就是m^2*s,
如果是弯矩功率谱密度，单位就是(N*m)^2*s
位移功率谱——m^2*s
速度功率谱——m^2/s
加速度功率谱——m^2/s^3

信号傅立叶变换的幅度图和频谱图matlab示例
fs=100;
x=-2:1/fs:2;
y=sin(3*pi*x);
z=rectpuls(x);
figure;plot(x,y,x,z,':r');

my=abs(fft(y));
mz=abs(fft(z));
my=my/max(my); %归一化
mz=mz/max(mz); %归一化

f=(0:1/length(x):1)*fs;
figure;plot(f(1:fs/2),my(1:fs/2),f(1:fs/2),mz(1:fs/2),':r');

my=20*log10(my+eps);mz=20*log10(mz+eps);
figure;plot(f(1:fs/2),my(1:fs/2),f(1:fs/2),mz(1:fs/2),':r'

时域信号--->相关函数--(FFT变换)-->功率谱--(除以频率分辨率)-->功率谱密度，这叫做间接求法，可以抑制白噪声，或者通俗的说不规律信号，分析的点数越多，规律信号的信噪比越好。

时域信号--(FFT变换)-->幅度谱--(平方)-->功率谱，这叫直接求法，最好不要用，除非你就想分析噪声有多

 

功率谱密度

　　在物理学中，信号通常是波的形式，例如电磁波、随机振动或者声波。当波的频谱密度乘以一个适当的系数后将得到每单位频率波携带的功率，这被称为信号的功率谱密度（power spectral density, PSD）或者谱功率分布（spectral power distribution, SPD）。功率谱密度的单位通常用每赫兹的瓦特数（W/Hz）表示，或者使用波长而不是频率，即每纳米的瓦特数（W/nm）来表示。
　　尽管并非一定要为信号或者它的变量赋予一定的物理量纲，下面的讨论中假设信号在时域内变化。
　　上面能量谱密度的定义要求信号的傅里叶变换必须存在，也就是说信号平方可积或者平方可加。一个经常更加有用的替换表示是功率谱密度（PSD），它定义了信号或者时间序列的功率如何随频率分布。这里功率可能是实际物理上的功率，或者更经常便于表示抽象的信号被定义为信号数值的平方，也就是当信号的负载为1欧姆(ohm)时的实际功率。此瞬时功率（平均功率的中间值）可表示为：
　　由于平均值不为零的信号不是平方可积的，所以在这种情况下就没有傅里叶变换。幸运的是维纳-辛钦定理（Wiener-Khinchin theorem）提供了一个简单的替换方法，如果信号可以看作是平稳随机过程，那么功率谱密度就是信号自相关函数的傅里叶变换。
　　信号的功率谱密度当且仅当信号是广义的平稳过程的时候才存在。如果信号不是平稳过程，那么自相关函数一定是两个变量的函数，这样就不存在功率谱密度，但是可以使用类似的技术估计时变谱密度。
　　f(t) 的谱密度和 f(t) 的自相关组成一个傅里叶变换对（对于功率谱密度和能量谱密度来说，使用着不同的自相关函数定义）。
　　通常使用傅里叶变换技术估计谱密度，但是也可以使用如Welch法（Welch's method）和最大熵这样的技术。
　　傅里叶分析的结果之一就是Parseval定理（Parseval's theorem），这个定理表明能量谱密度曲线下的面积等于信号幅度平方下的面积，总的能量是：
　　：上面的定理在离散情况下也是成立的。另外的一个结论是功率谱密度下总的功率与对应的总的平均信号功率相等，它是逐渐趋近于零的自相关函数。
　　功率谱密度谱是一种概率统计方法，是对随机变量均方值的量度。一般用于随机振动分析，连续瞬态响应只能通过概率分布函数进行描述，即出现某水平响应所对应的概率。
　　功率谱密度的定义是单位频带内的“功率”（均方值）
　　功率谱密度是结构在随机动态载荷激励下响应的统计结果，是一条功率谱密度值—频率值的关系曲线，其中功率谱密度可以是位移功率谱密度、速度功率谱密度、加速度功率谱密度、力功率谱密度等形式。数学上，功率谱密度值—频率值的关系曲线下的面积就是方差，即响应标准偏差的平方值。

 

 

功率谱密度

在物理学中，信号通常是波的形式，例如电磁波、随机振动或者声波。当波的频谱密度乘以一个适当的系数后将得到每单位频率波携带的功率，这被称为信号的功率谱密度（power spectral density, PSD）或者谱功率分布（spectral power distribution, SPD）。功率谱密度的单位通常用每赫兹的瓦特数（W/Hz）表示，或者使用波长而不是频率，即每纳米的瓦特数（W/nm）来表示。