您好,欢迎访问三七文档
t=0:0.0001:0.1;%时间间隔为0.0001,说明采样频率为10000Hzx=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采样点的个数为1000subplot(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));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]);改进的直接法:对于直接法的功率谱估计,当数据长度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);做FFT分析时,幅值大小与FFT选择的点数有关,但不影响分析结果。在IFFT时已经做了处理。要得到真实的振幅值的大小,只要将得到的变换后结果乘以2除以N即可。自相关x(n)和y(n)为统计的随机序列,-∞n∞,其中E{·}为预期的数值操作,且xcorr函数只能估计有限序列。c=xcorr(x,y)返回矢量长度为2*N-1互相关函数序列,其中x和y的矢量长度均为N,如果x和y的长度不一样,则在短的序列后补零直到两者长度相等输出矢量C通过c(m)=Rxy(m-N),m=1,...,2N-1.得到。通常,互相关函数需要正规化来得到更准确的估计值。c=xcorr(x)为矢量x的自相关估计;c=xcorr(x,y,'option')为有正规化选项的互相关计算;其中选项为biased为有偏的互相关函数估计;unbiased为无偏的互相关函数估计;coeff为0延时的正规化序列的自相关计算;none为原始的互相关计算;若需了解更多的有偏和无偏相关估计,请参阅文献[1].c=xcorr(x,'option')特指以上某个选项的自相关估计。c=xcorr(x,y,maxlags)返回一个延迟范围在[-maxlags,maxlags]的互相关函数序列,输出c的程度为2*maxlags+1.c=xcorr(x,maxlags)返回一个延迟范围在[-maxlags,maxlags]的自相关函数序列,输出c的程度为2*maxlags+1.c=xcorr(x,y,maxlags,'option')同时指定maxlags和option的互相关计算.c=xcorr(x,maxlags,'option')同时指定maxlags和option的自相关计算.[c,lags]=xcorr(...)返回一个在c进行相关估计的延迟矢量lag,其范围为[-maxlags:maxlags],当maxlags没有指定时,其范围为[-N+1,N-1]lags就是信号的延时或者超前啊,两个信号的相关性是相对于一定超前和滞后而言的。相关算法就是用移位相乘来体现信号之间的相似度,这里包括幅值和频率。你看看帮助里面的例子:x=[1,2i,3];y=[4,5,6];[c1,lags]=xcorr(x,y)
本文标题:功率谱密度
链接地址:https://www.777doc.com/doc-2611979 .html