1、求信号功率谱时候用下面的不同方法，功率谱密度的幅值大小相差很大！我的问题是，计算具体信号时，到底应该以什么准则决定该选用什么方法啊？功率谱密度的幅植的具体意义是什么？下面是一些不同方法计算同一信号的 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*
2、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)+
3、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)分
4、段求周期图再平均。 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_Px
5、xc=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*p
6、i*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=pw
7、elch(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);-功率谱的数据都是相对值，他无法给出信号的实际绝对幅值，一般只要看峰值之间的比值正确就行了，当然这个问题可以通过做正规化处理解决-谢谢回答!不过具体原因,我还是不很明白
8、啊?你能不能从原理上讲讲看啊! 而且有时候所求信号的幅值意义是很大的!比如,地震信号的功率谱值,其幅值有一定的范围,而我求出来的值总是和文献的对不上,不知道具体选择求解方法时怎么处理啊 ?-xn=cos(2*pi*40*n)+3*cos(2*pi*100*n)像这个信号的话,考虑单边谱的话,40Hz 的幅值为 1,100Hz 的幅值为 3,对应的功率谱值分别为 1 和 9.在用 FFT 做谱估计值时,应把 FFT 的结果取模后除以 FFT 的点数再乘以 2,得到单边谱幅值,再平方后就得到单边功率谱值-大家继续讨论啊!在地震信号处理中，常常统计功率谱幅值的变化规律，按大家的说法，就毫无意义了，因
9、为采用不同方法，幅值差别很大，到底是怎么一回事啊？-实际上，功率谱也可以求出真实幅值的，只要把求得的结果对采样频率归一化，即按楼主的方法得到的相对值再乘以采样频率即为真实结果！-用各种方法得到的功率谱，都是相对值，往往用分贝来表示。则是否能知道它的绝对值，这是可能的，但需要进行校正，要设定 0 分贝对应的数值。我们测量的无非是加速度、速度和位移的功率谱，则它们 0 分贝的参考值分别是：a0=1um/s2v0=1nm/sd0=1pm（摘自 马大猷 声学手册）。为了能求出绝对的功率谱密度数值，要从测量源头开始便要进行校正，也就是从传感器、放大器、直至进 AD 变换，每一环节都要校正，这样才能拿到正
10、确的结果。-功率谱密度，单位为：unit2/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 乘以
11、0.5，正好为该信号有效值 0.707 的平方。-能不能解释一下为什么要转化为单边谱，而且最后除 0.5 是什么意思？谢谢-因为实数信号的双边幅值谱都是对称的，因此用单边谱就够了，这时候把负频率成分附加到相应的正频率成分，也就是在双边谱的基础上乘以双边谱中包含负频率，在物理系统中是没有的；0.5 为频率分辨率。-引用:原帖由 yangzp 于 2006-9-23 15:31 发表功率谱密度，单位为：unit2/Hz 代表单位频率上信号的能量，所以是密度谱，幅值代表频段内的有效值平方，计算时的步骤为1 对每一分段数据进行 FFT 变换，并求的幅值谱2 对幅值谱进行平方3 将双边谱转化为单边谱 .
12、 我周期图法得到的幅值为 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
13、索引+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
14、+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 的结果取
15、模后除以 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 的平方。大家看，他们两个说的明显的差
16、一个二倍，到底那个是对的呢？-1)求信号功率谱时候用下面的不同方法，功率谱密度的幅值大小相差很大！matlab 提供的功率谱求法肯定是没有错的，关键在于各个参数的选取，以及对每种求解方法的了解！平均周期图法和其他方法求出的结果，参数条件取得一样的话，可以得到完全相同的结果。2）我的问题是，计算具体信号时，到底应该以什么准则决定该选用什么方法啊？直接采用平均周期图法求功率谱时，功率普形状呈锯齿形，谱峰点的准确位置不大好定。于是可以采用其他的方法对谱进行平滑操作，平滑化，仅仅是为了使图形光滑，并不会使得波的本质受到歪曲和畸变。反过来说，由于不纯的东西去掉了，本质的东西必然会更加显示出来！平滑化的程
17、度可以根据所分析的信号，选择合理的窗函数和带宽！可以采用逐步观察的方法进行考察！ 3）功率谱密度的幅植的具体意义是什么？对于地震动信号（我研究的范畴），对于功率谱的单位，如数据为加速度时，为米（平方）/ 秒（三次方）；若数据为速度或者位移时，他的单位为米（平方）/秒，米（平方）*秒；他和物理学中的功率（单位时间内所作的功，单位为瓦特）概念是不一样的。因此在地震波的情况下，所谓功率谱，不能用功率这一物理量来理解。对于一个问题的理解，关键靠自己多实践，找相关的资料学习，这样可以加深理解-我现在通过试验得到离散化的峰值谱曲线，按上面讨论的提示，是不是只要先进行平方，然后除以频率分辨率就可以了？我们试
18、验测试的都是单边谱，是否还需要2。谢谢大家解答。-你肯定理解错了单边谱和双边谱的概念了。单边谱就是双边谱乘以 2，所以你得到的如果是单边谱的话，那么应该在平方然后除以频率分辨率后，再除以 2-引用 : 原帖由 yangzp 于 2006-9-23 15:31 发表 功率谱密度，单位为：unit2/Hz 代表单位频率上信号的能量，所以是密度谱，幅值代表频段内的有效值平方，计算时的步骤为1 对每一分段数据进行 FFT 变换，并求的幅值谱2 对幅值谱进行平方3 将双边谱转化为单边谱. 按Pxx,f=pwelch(xn,window,noverlap,nfft,fs,range);求得的 Pxx 是不
19、是等于 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
20、 点 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(s
21、ize(n); nfft=1024; Pxx=fft(xn).*conj(fft(xn);%直接法 plot(10*log10(Pxx),b); holdFs=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);图中 直接法功率谱(兰) 间接法功率谱(红)
