功率谱密度
- 1、下载文档前请自行甄别文档内容的完整性,平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益,请联系客服反馈,我们会尽快为您处理(人工客服工作时间:9:00-18:30)。
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点数据的傅里叶变换(即频谱),然后取频谱与其共轭的乘积,就得到功率谱的估计;间接法先计
算N点样本数据的自相关函数,然后取自相关函数的傅里叶变换,即得到功率谱的估计.都可以编程实现,很简单。在matlab中,周期图法可以用函数periodogram实现。
但是周期图法估计出的功率谱不够精细,分辨率比较低。因此需要对周期图法进行修正,可以将信号序列x(n)分为n个不相重叠的小段,分别用周期图法进行谱估计,然后将这n段数据估计的结果的平均值作为整段数据功率谱估计的结果。还可以将信号序列x(n)重叠分段,分别计算功率谱,再计算平均值作为整段数据的功率谱估计。
这2种称为分段平均周期图法,一般后者比前者效果好。加窗平均周期图法是对分段平均周期图法的改进,即在数据分段后,对每段数据加一个非矩形窗进行预处理,然后在按分段平均周期图法估计功率谱。相对于分段平均周期图法,加窗平均周期图法可以减小频率泄漏,增加频峰的宽度。welch法就是利用改进的平均周期图法估计估计随机信号的功率谱,它采用信号分段重叠,加窗,FFT等技术来计算功率谱。与周期图法比较,welch法可以改善估计谱曲线的光滑性,大大提高谱估计的分辨率。
现代谱估计主要针对经典谱估计分辨率低和方差性不好提出的,可以极大的提高估计的分辨率和平滑性。可以分为参数模型谱估计和非参数模型谱估计。参数模型谱估计有AR模型,MA模型,ARMA模型等;非参数模型谱估计有最小方差法和MUSIC法等。由于涉及的问题太多,这里不再详述,可以参考有关资料。
matlab中,现代谱估计的很多方法都可以实现。music方法用pmusic 命令实现;pburg函数利用burg法实现功率谱估计;pyulear函数利用yule-walker算法实现功率谱估计等等。
另外,sptool工具箱也具有功率谱估计的功能。窗口化的操作界面很方便,而且有多种方法可以选择。
经典功率谱估计
直接法:
直接法又称周期图法,它是把随机序列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));