matlab 功率谱分析

1、直接法：
: }' |$ |4 Y* o, H直接法又称周期图法，它是把随机序列x(n)的N个观测数据视为一能量有限的序列，直接计算x(n)的离散傅立叶变换，得X(k)，然后再取其幅值的平方，并除以N，作为序列x(n)真实功率谱的估计。; ?# C. [, B5 ^( S# M: P- X
Matlab代码示例：
9 d" ^! M0 c  r7 J  p: D. pclear;: p/ k; x  W% O
Fs=1000; %采样频率
4 _3 C6 K2 N( F4 n2 Yn=0:1/Fs:1;3 H+ y* Z/ ~& Y- v
%产生含有噪声的序列
0 Z0 n1 s/ C0 Z* H3 Axn=cos(2*pi*40*n)+3*cos(2*pi*100*n)+randn(size(n));
8 T& k$ E! n* [8 Cwindow=boxcar(length(xn)); %矩形窗
) i% Y$ f$ B/ M) V" i5 }nfft=1024;
4 s" L% [0 T/ Y; X$ G4 E[Pxx,f]=periodogram(xn,window,nfft,Fs); %直接法
/ d! _! g0 N  [8 ?) S6 K! Cplot(f,10*log10(Pxx));: O8 M& X; @8 H+ Z& O" _
; N' g; O+ R, s
2、间接法：( a  a- o8 K; F; L$ j# U: k! g
间接法先由序列x(n)估计出自相关函数R(n)，然后对R(n)进行傅立叶变换，便得到x(n)的功率谱估计。
' D/ l) i+ O2 `9 EMatlab代码示例：( a+ j. p- S* B# [+ K& b
clear;
* m2 ?* ?' R3 sFs=1000; %采样频率
) f; a# A* k( e8 [2 A& b; ?9 Zn=0:1/Fs:1;
( P. s, Q6 I( C# N1 N%产生含有噪声的序列, B! ]) y2 F' l, B
xn=cos(2*pi*40*n)+3*cos(2*pi*100*n)+randn(size(n));" C. I* B7 B. u- f
nfft=1024;
  T4 u2 u# a& z9 l- [7 |8 ?cxn=xcorr(xn,'unbiased'); %计算序列的自相关函数
- K/ @" ?% h& I& A" ^' {, b1 \1 `CXk=fft(cxn,nfft);" R" B! G8 L- c4 W7 t2 b
Pxx=abs(CXk);4 k9 g0 `- k: @, C+ U8 I  q2 M
index=0:round(nfft/2-1);
: O5 s  t& S9 q* l! G+ m. i6 X) Bk=index*Fs/nfft;
* t( v% _% V& S$ ^. dplot_Pxx=10*log10(Pxx(index+1));
7 @! L- \& r7 Tplot(k,plot_Pxx);( k- u' W( ~; Z) D+ F" t% v! A; W* y

" @/ c/ V, u6 K: H* I* q+ o" u3、改进的直接法：) k: b: r7 S8 H
对于直接法的功率谱估计，当数据长度N太大时，谱曲线起伏加剧，若N太小，谱的分辨率又不好，因此需要改进。+ k  a2 E3 q0 }% z+ d$ @
3.1、Bartlett法
5 O$ `. R  Z# B2 DBartlett平均周期图的方法是将N点的有限长序列x(n)分段求周期图再平均。- L- f8 U0 N3 j. j
Matlab代码示例：
: L6 a! K7 }+ A! ~, p% sclear；
. d6 I- O7 u) }/ Q; W& IFs=1000;  F2 a& c( P- g, _5 I7 W
n=0:1/Fs:1;
8 d, I" r' Z7 }" `xn=cos(2*pi*40*n)+3*cos(2*pi*100*n)+randn(size(n));; U, R+ p: T2 V4 ?7 S  P, h# U5 }
nfft=1024;* d, c* O% O5 |2 B9 l# d! L) B
window=boxcar(length(n)); %矩形窗' w+ O6 o" m9 d/ I8 o3 `% @
noverlap=0; %数据无重叠
* i$ x5 w( {" P/ i2 A+ |p=0.9; %置信概率( x5 T$ _1 u; H: q0 d8 q
[Pxx,Pxxc]=psd(xn,nfft,Fs,window,noverlap,p);
0 S, a5 a/ M' I0 d% W% N" z" Kindex=0:round(nfft/2-1);
' b% G% f7 M  S. n* N- J5 L8 bk=index*Fs/nfft;9 Y4 _0 \) V4 K. l" T$ u
plot_Pxx=10*log10(Pxx(index+1));
0 l% F9 B8 _0 H9 K) u: z6 u) yplot_Pxxc=10*log10(Pxxc(index+1));
7 I$ f6 h( g$ v0 K# ^& Ifigure(1)
( a, c; V! O3 G: w5 ?plot(k,plot_Pxx);0 C6 P1 K, ]+ @
pause;8 i9 _1 s: a2 Y( A& |, q) i: ^
figure(2)
3 ~4 W! e" L7 g  splot(k,[plot_Pxx plot_Pxx-plot_Pxxc plot_Pxx+plot_Pxxc]);
* b# G/ q" a+ V3 N1 V) |% h
9 e6 [  g9 e7 e- ]' j# B3.2、Welch法
0 P9 A! F1 Q+ s/ D/ M! mWelch法对Bartlett法进行了两方面的修正，一是选择适当的窗函数w(n)，并再周期图计算前直接加进去，加窗的优点是无论什么样的窗函数均可使谱估计非负。二是在分段时，可使各段之间有重叠，这样会使方差减小。 $ @4 O4 @  f& u9 {8 p& I
Matlab代码示例：, W! o9 b" s* u7 c# i
clear;
" b2 f6 Z  `! |$ h0 }Fs=1000;7 \4 H4 I& {0 M
n=0:1/Fs:1;. ^* }: n, f% x4 s" l: E2 q2 r
xn=cos(2*pi*40*n)+3*cos(2*pi*100*n)+randn(size(n));
( z0 C7 t9 r+ m$ f- Bnfft=1024;) i* ~$ e" X" S. {, ^
window=boxcar(100); %矩形窗
6 [0 W3 G+ e& r! R7 y% jwindow1=hamming(100); %海明窗
0 o2 E& _& h8 Y* n; T$ I: Jwindow2=blackman(100); %blackman窗
1 Y0 t; w' D. Y; Dnoverlap=20; %数据无重叠
7 \5 i7 i2 v+ ^6 ]+ Srange='half'; %频率间隔为[0 Fs/2]，只计算一半的频率
2 q! y" b- J( O[Pxx,f]=pwelch(xn,window,noverlap,nfft,Fs,range);
8 d+ f  W' s/ @5 l[Pxx1,f]=pwelch(xn,window1,noverlap,nfft,Fs,range);
  S. z3 u! X* T' x- F+ A8 @# q[Pxx2,f]=pwelch(xn,window2,noverlap,nfft,Fs,range);
! a; ~1 r1 @8 c2 v8 J0 `; [plot_Pxx=10*log10(Pxx);
: U; e' ]7 ~9 Y- q) H  p8 {plot_Pxx1=10*log10(Pxx1);
; R( ~+ i* I' F, Oplot_Pxx2=10*log10(Pxx2);" Q- H3 g$ P4 l) x2 D6 d! C& Y; x

  {) |; z; s8 \figure(1)4 g6 [+ q5 b% I+ ]7 i
plot(f,plot_Pxx);. X3 O: a: @) G& A" Z

1 Q7 h- T* A- t, B0 }; Wpause;
0 N* h0 N3 t' Y7 ?8 R6 ~9 ~0 n/ g" A
figure(2)) [3 o0 j/ s) S7 r4 F
plot(f,plot_Pxx1);
' ?( z9 u: U5 a& {
( E0 E$ K. [. |# ]& N! rpause;
8 n. z" g, q" ^' j
. g2 o) l7 y1 D+ \1 I- t/ S0 ?figure(3)6 r1 Q; R1 s2 f  W6 A
plot(f,plot_Pxx2);