matlab 功率谱分析

1、直接法：
. l- L# c8 [* l4 ?4 \' ?直接法又称周期图法，它是把随机序列x(n)的N个观测数据视为一能量有限的序列，直接计算x(n)的离散傅立叶变换，得X(k)，然后再取其幅值的平方，并除以N，作为序列x(n)真实功率谱的估计。# R  d; l5 O  ]7 ~( w9 c4 |' e
Matlab代码示例：( b* i1 i- p' ?2 j
clear;8 P# Y+ V& |$ Y$ q& {3 }
Fs=1000; %采样频率1 H$ B  v/ |+ j8 ^7 t1 Z$ k
n=0:1/Fs:1;
; p* s- {! `  [$ d4 Q8 x; S/ ~, C) Y%产生含有噪声的序列: Q* R! @* m% R' o7 z
xn=cos(2*pi*40*n)+3*cos(2*pi*100*n)+randn(size(n));
% x# |2 I2 d& v2 x3 s% S; h' Jwindow=boxcar(length(xn)); %矩形窗, h. Z) a# Z; m* v* C0 w* k
nfft=1024;
; o! D5 }" O: p; v( I[Pxx,f]=periodogram(xn,window,nfft,Fs); %直接法" v. g( H: _8 _" W( ?8 y
plot(f,10*log10(Pxx));# s( p/ }8 Q  g% Y: d* O5 {

- {5 g* z- z; \+ V6 |$ j2、间接法：
3 n- @5 ^! i* p. e7 p4 h' h间接法先由序列x(n)估计出自相关函数R(n)，然后对R(n)进行傅立叶变换，便得到x(n)的功率谱估计。
2 I# \! ?5 |! |, kMatlab代码示例：5 F. I/ O6 R9 ^% F6 {3 u
clear;" i8 w) {) z" N1 v
Fs=1000; %采样频率$ J9 I0 K9 t; B% |$ [4 x2 z
n=0:1/Fs:1;$ {, V. P% [2 Z0 I" n4 y
%产生含有噪声的序列' B: t2 ]1 y9 v, o5 m% s! i
xn=cos(2*pi*40*n)+3*cos(2*pi*100*n)+randn(size(n));
5 s8 R0 X! S: P, G( x& {nfft=1024;
4 @8 k% V4 T6 P& `6 ?- G$ hcxn=xcorr(xn,'unbiased'); %计算序列的自相关函数
. P" ~8 u; ]7 Q9 @& |4 aCXk=fft(cxn,nfft);
6 E, B. x3 P2 A6 o1 f+ XPxx=abs(CXk);: @9 s" o1 {& u5 p0 O0 _3 F: b5 ~
index=0:round(nfft/2-1);' u. w: ~7 @/ y6 i" d+ i! }' N
k=index*Fs/nfft;
5 _$ Z) ?, d3 }5 g3 X: e7 Kplot_Pxx=10*log10(Pxx(index+1));
8 X* E8 ]! {: h6 k) `3 m, K0 D! uplot(k,plot_Pxx);, q1 C5 c! r5 f

/ M1 h8 B( i: b+ A4 P3、改进的直接法：
  s  P( l( K/ [8 D6 h对于直接法的功率谱估计，当数据长度N太大时，谱曲线起伏加剧，若N太小，谱的分辨率又不好，因此需要改进。
6 q2 ~3 P1 D: Q3.1、Bartlett法5 `# D  x& _0 c1 \" v
Bartlett平均周期图的方法是将N点的有限长序列x(n)分段求周期图再平均。
+ i" M4 t% f. H+ ~Matlab代码示例：0 S8 ^3 p- f. W
clear；
* V+ {/ S- y" X. R( OFs=1000;
1 j; r; q8 g' t! U9 \5 N& zn=0:1/Fs:1;& L- ~( k: w) R5 e$ C# b; D3 O
xn=cos(2*pi*40*n)+3*cos(2*pi*100*n)+randn(size(n));/ l( V. c: K2 a. i& b! Q& U9 c
nfft=1024;5 f# G! s' ?- e$ @- r
window=boxcar(length(n)); %矩形窗) v* n3 Z6 d! L- q  V& K& z
noverlap=0; %数据无重叠
% L8 F/ R5 g  p6 E5 W2 Mp=0.9; %置信概率4 D8 h& Y- V( Z2 I+ @2 m
[Pxx,Pxxc]=psd(xn,nfft,Fs,window,noverlap,p);- U2 e/ A' f  e4 w1 z0 P* f( ^
index=0:round(nfft/2-1);8 G4 T2 ]8 B0 f, X4 A
k=index*Fs/nfft;
6 x- E! _9 h6 z1 |4 F' ]plot_Pxx=10*log10(Pxx(index+1));
& U5 U! H9 v! Y5 W1 a; g/ Qplot_Pxxc=10*log10(Pxxc(index+1));
* r+ [7 P, @4 Dfigure(1)/ Q4 m$ u$ Y6 o. \( \$ M3 u
plot(k,plot_Pxx);8 i4 S) H" P6 `. U+ _
pause;
; D. E' k3 S' j% yfigure(2)
! n  Q5 ?* v9 _  Q* x- W! }+ Qplot(k,[plot_Pxx plot_Pxx-plot_Pxxc plot_Pxx+plot_Pxxc]);6 x: T1 T: q& q6 F; I
4 D$ X  j- J0 o; Q2 p+ @' F6 j
3.2、Welch法
% s9 [- U) ~& [Welch法对Bartlett法进行了两方面的修正，一是选择适当的窗函数w(n)，并再周期图计算前直接加进去，加窗的优点是无论什么样的窗函数均可使谱估计非负。二是在分段时，可使各段之间有重叠，这样会使方差减小。 + b4 R' r; g. q$ u, }9 Y& l
Matlab代码示例：
' n8 q5 Q) N: z7 jclear;. T1 g" Q7 ^. j* J' x3 N
Fs=1000;, |5 h! K2 c6 o3 V& _
n=0:1/Fs:1;& e: T" r/ |0 e
xn=cos(2*pi*40*n)+3*cos(2*pi*100*n)+randn(size(n));
5 L9 b$ z* ?  j7 X; @nfft=1024;. f# C$ k3 W+ u
window=boxcar(100); %矩形窗4 X- `# R  M1 ]
window1=hamming(100); %海明窗
, _$ i! @* B0 d1 S- Rwindow2=blackman(100); %blackman窗! y, X  p  s! z3 [0 _
noverlap=20; %数据无重叠
+ R" h' j1 a0 [+ M& m3 j# z! n0 Wrange='half'; %频率间隔为[0 Fs/2]，只计算一半的频率
6 w9 |5 t6 L. P6 s' a[Pxx,f]=pwelch(xn,window,noverlap,nfft,Fs,range);
  W0 d* d5 o2 n+ a) c! x[Pxx1,f]=pwelch(xn,window1,noverlap,nfft,Fs,range);
6 p( |& P8 n3 Y9 Z6 ~5 P8 V[Pxx2,f]=pwelch(xn,window2,noverlap,nfft,Fs,range);) x/ p5 u5 }( N2 z. O
plot_Pxx=10*log10(Pxx);
6 Q8 C" _- v0 S5 K& v' n. splot_Pxx1=10*log10(Pxx1);& H, f( V0 b3 v/ J
plot_Pxx2=10*log10(Pxx2);
" e( f; @8 B, W$ B7 D8 }/ b4 l; [
figure(1)
, a8 o) e, ?; ^- Lplot(f,plot_Pxx);& T7 e' M' [, s. \/ {" |) r
. }5 I+ `5 ^& Y" m6 Q# q
pause;% u5 o% `1 u' B
+ Q5 C2 G5 Z9 W% f5 r* o' o7 c
figure(2)
6 G! n) `8 O! P: U% g, A" h( R( wplot(f,plot_Pxx1);
( s- v+ y. q4 R' j- f! @) h1 j& O& n% x' X) z1 N
pause;
: K) A4 ~' S; p2 j+ W! d& J
; i- _1 b$ B2 _* [' |figure(3)
" ~# g  l# B, r& ]plot(f,plot_Pxx2);