matlab 功率谱分析

mutougeda

matlab 功率谱分析

1、直接法：
直接法又称周期图法，它是把随机序列x(n)的N个观测数据视为一能量有限的序列，直接计算x(n)的离散傅立叶变换，得X(k)，然后再取其幅值的平方，并除以N，作为序列x(n)真实功率谱的估计。
! `6 Y- j; L3 v( p9 CMatlab代码示例：
clear;
: u2 J! m! L$ B7 m8 r& @Fs=1000; %采样频率
, F/ Y* C5 I8 a1 p2 _# v; O! dn=0:1/Fs:1;
%产生含有噪声的序列
xn=cos(2*pi*40*n)+3*cos(2*pi*100*n)+randn(size(n));
/ S) ^% v+ Y7 ~6 D3 _8 Mwindow=boxcar(length(xn)); %矩形窗
9 g  e3 ?4 E% V+ \  e$ qnfft=1024;
4 v9 g- f' |+ V. j% c[Pxx,f]=periodogram(xn,window,nfft,Fs); %直接法
+ z* a6 b4 t7 Z/ o2 k& P# splot(f,10*log10(Pxx));

2、间接法：
/ C! X% A# l: i) F$ U# s& `3 x0 D间接法先由序列x(n)估计出自相关函数R(n)，然后对R(n)进行傅立叶变换，便得到x(n)的功率谱估计。
1 G3 s, ?2 c  V& k  i4 ?# C$ F1 PMatlab代码示例：
clear;
  A+ b" D5 {: G  h- X8 @1 h) v" GFs=1000; %采样频率
; m5 b, {& y2 \( e/ P3 dn=0:1/Fs:1;
" ~- {6 m- [- B. p8 x: T# L%产生含有噪声的序列
xn=cos(2*pi*40*n)+3*cos(2*pi*100*n)+randn(size(n));
" I, ^7 O& T, j0 S) Jnfft=1024;
8 ?* T9 P! M- j  hcxn=xcorr(xn,'unbiased'); %计算序列的自相关函数
9 |" z% ?7 W9 q. `( SCXk=fft(cxn,nfft);
Pxx=abs(CXk);
index=0:round(nfft/2-1);
k=index*Fs/nfft;
5 \8 W  ^4 ?0 F# L7 \plot_Pxx=10*log10(Pxx(index+1));
- l4 A. B2 L8 v2 Jplot(k,plot_Pxx);

0 g+ E5 ~1 Z3 t5 d& q% s3、改进的直接法：
对于直接法的功率谱估计，当数据长度N太大时，谱曲线起伏加剧，若N太小，谱的分辨率又不好，因此需要改进。
1 F9 Q% Q& K6 h* `$ ]+ J) o; [3.1、Bartlett法
/ T' t2 H' g" b% Q7 L$ QBartlett平均周期图的方法是将N点的有限长序列x(n)分段求周期图再平均。
! I, G  m  I3 q6 OMatlab代码示例：
clear；
Fs=1000;
% s0 }3 `4 C( T8 C5 E1 J6 s. Ln=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);
! l8 A0 C( |+ J* Y7 {k=index*Fs/nfft;
) q$ o/ e& F, S7 m" oplot_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]);

3.2、Welch法
/ a/ M! P! F% R. F/ O4 @. qWelch法对Bartlett法进行了两方面的修正，一是选择适当的窗函数w(n)，并再周期图计算前直接加进去，加窗的优点是无论什么样的窗函数均可使谱估计非负。二是在分段时，可使各段之间有重叠，这样会使方差减小。
, |4 G; `, i: }% j5 q' j7 SMatlab代码示例：
8 T) M% `# v, F* t0 u$ l, H, p) Xclear;
  |% |. z: p" z( n* K  O, W$ hFs=1000;
6 k7 X! S8 B6 En=0:1/Fs:1;
) G' f: a6 I& t8 J; kxn=cos(2*pi*40*n)+3*cos(2*pi*100*n)+randn(size(n));
3 H7 i$ n. X4 ~# r8 k$ B/ Dnfft=1024;
+ h3 i! _6 X( fwindow=boxcar(100); %矩形窗
8 \  H5 Z* N# |5 T" z2 `6 H1 Xwindow1=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);

' G$ |3 i+ @0 ]- Z) Mfigure(1)
8 ~. S( @( r3 O4 b- Kplot(f,plot_Pxx);
8 h, L" k2 C4 R9 P. d
6 o8 i- z7 @) opause;
% E; w; O4 W  R
figure(2)
/ B5 Y, I9 m3 Y$ m4 G% Fplot(f,plot_Pxx1);

pause;

8 @5 g3 P9 {$ t: d6 w& y- Wfigure(3)
+ S+ x3 _7 t; o$ Qplot(f,plot_Pxx2);

2 N$ q9 N8 l, b4 W- i& l1 \3 |

NingW

matlab 功率谱分析

Heaven_1

plot函数比较好