求信号功率谱时候用下面的不同方法.docx
《求信号功率谱时候用下面的不同方法.docx》由会员分享,可在线阅读,更多相关《求信号功率谱时候用下面的不同方法.docx(9页珍藏版)》请在冰豆网上搜索。
求信号功率谱时候用下面的不同方法
求信号功率谱时候用下面的不同方法,功率谱密度的幅值大小相差很大!
我的问题是,计算具体信号时,到底应该以什么准则决定该选用什么方法啊?
功率谱密度的幅植的具体意义是什么?
?
下面是一些不同方法计算同一信号的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_Pxxplot_Pxx-plot_Pxxcplot_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';%频率间隔为[0Fs/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是什么意思?
谢谢
------------------------------
因为实数信号的双边幅值谱都是对称的,因此用单边谱就够了,这时候把负频率成分附加到相应的正频率成分,也就是在双边谱的基础上乘以2.
双边谱中包含负频率,在物理系统中是没有的;0.5为频率分辨率。
-----------------------------
引用:
原帖由yangzp于2006-9-2315:
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-2315:
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后计功率谱,
closeall;clc;clearall;
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)+[0cxn(1:
nfft-1)];
CXk=fft(cxn,nfft);
Pxx=abs(CXk);
plot_Pxx=10*log10(Pxx);
plot(plot_Pxx,'r');
图中直接法功率谱(兰) 间接法功率谱(红)