求信号功率谱时候用下面不同方法

合集下载
相关主题
  1. 1、下载文档前请自行甄别文档内容的完整性,平台不提供额外的编辑、内容补充、找答案等附加服务。
  2. 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
  3. 3、如文档侵犯您的权益,请联系客服反馈,我们会尽快为您处理(人工客服工作时间:9:00-18:30)。

求信号功率谱时候用下面的不同方法,功率谱密度的幅值大小相差很大!

我的问题是,计算具体信号时,到底应该以什么准则决定该选用什么方法啊?

功率谱密度的幅植的具体意义是什么??下面是一些不同方法计算同一信号的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=; %置信概率

[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,得到单边谱幅值,再平方后

就得到单边功率谱值

--------------------------------

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

相关文档
最新文档