matlab 功率谱分析

mutougeda

matlab 功率谱分析

1、直接法：
直接法又称周期图法，它是把随机序列x(n)的N个观测数据视为一能量有限的序列，直接计算x(n)的离散傅立叶变换，得X(k)，然后再取其幅值的平方，并除以N，作为序列x(n)真实功率谱的估计。
& d# q: @* @  ^. hMatlab代码示例：
clear;
, o5 F3 j2 T& e9 Z1 ~- OFs=1000; %采样频率
n=0:1/Fs:1;
%产生含有噪声的序列
: P9 F) k  s6 w5 n; J7 Ixn=cos(2*pi*40*n)+3*cos(2*pi*100*n)+randn(size(n));
% I0 E. [, M  _& [0 Wwindow=boxcar(length(xn)); %矩形窗
  h. t0 U& Y, K% c1 Q: b% y7 Ynfft=1024;
8 u( S' O# N* g8 d8 c6 j, G' {7 n# W# u% r[Pxx,f]=periodogram(xn,window,nfft,Fs); %直接法
( f- x8 i9 u; W8 X/ a: ~. C3 pplot(f,10*log10(Pxx));

  t2 {9 B/ g: q2、间接法：
间接法先由序列x(n)估计出自相关函数R(n)，然后对R(n)进行傅立叶变换，便得到x(n)的功率谱估计。
5 t8 V7 u& ]( F7 BMatlab代码示例：
2 x+ n4 b) E8 x2 ]$ Jclear;
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);
8 B* l/ W: d- S* M6 B1 A. }$ mPxx=abs(CXk);
index=0:round(nfft/2-1);
k=index*Fs/nfft;
0 ~. d; L4 O$ o# A: Qplot_Pxx=10*log10(Pxx(index+1));
" k0 G& z+ R" s. Mplot(k,plot_Pxx);
8 n2 p, x4 c7 t* l. y3 Y
/ Z* e; o6 l* n3、改进的直接法：
. L, [0 ?/ P, O6 K1 f6 d2 s对于直接法的功率谱估计，当数据长度N太大时，谱曲线起伏加剧，若N太小，谱的分辨率又不好，因此需要改进。
3.1、Bartlett法
Bartlett平均周期图的方法是将N点的有限长序列x(n)分段求周期图再平均。
Matlab代码示例：
# s( a& W* j: P) f! ?clear；
2 x3 j: w: A; E& H9 m% V" q2 |5 RFs=1000;
n=0:1/Fs:1;
$ c! [0 O/ B5 R. x5 M0 H$ Pxn=cos(2*pi*40*n)+3*cos(2*pi*100*n)+randn(size(n));
: T, [3 w: F* ]* ~5 ~$ M3 Ynfft=1024;
* [) d6 P3 e' Q5 vwindow=boxcar(length(n)); %矩形窗
; y, Z3 _6 H, jnoverlap=0; %数据无重叠
# b; M1 S# V  p$ C. E0 P$ ]p=0.9; %置信概率
0 g  r4 @; R+ N* O, I[Pxx,Pxxc]=psd(xn,nfft,Fs,window,noverlap,p);
7 ~) J. _& i. Z& F( findex=0:round(nfft/2-1);
k=index*Fs/nfft;
, u! H3 R. A2 ?1 {( Nplot_Pxx=10*log10(Pxx(index+1));
plot_Pxxc=10*log10(Pxxc(index+1));
0 h3 n; w+ B6 v. {# Ifigure(1)
  |2 ?" i7 |5 O: {6 z, N) Y- P, Q- M0 uplot(k,plot_Pxx);
pause;
: ~: @/ `  i- wfigure(2)
8 o7 v0 x! c) P  ~1 Z7 |/ }3 c# splot(k,[plot_Pxx plot_Pxx-plot_Pxxc plot_Pxx+plot_Pxxc]);
3 K, I3 a: a+ O- [/ l6 a$ F
! j! c$ U6 Z) e9 E% z: P3 V3.2、Welch法
/ L& S2 |" Z& h% ]; ZWelch法对Bartlett法进行了两方面的修正，一是选择适当的窗函数w(n)，并再周期图计算前直接加进去，加窗的优点是无论什么样的窗函数均可使谱估计非负。二是在分段时，可使各段之间有重叠，这样会使方差减小。
Matlab代码示例：
. p7 A( k' Z$ O2 B- lclear;
9 l) V/ _, b( HFs=1000;
2 w7 u0 O  x% ^' _8 H: gn=0:1/Fs:1;
. o2 r+ |6 ^# {6 Uxn=cos(2*pi*40*n)+3*cos(2*pi*100*n)+randn(size(n));
nfft=1024;
( x9 f* G4 z- a1 uwindow=boxcar(100); %矩形窗
& b1 ]3 P0 ]) u0 ?; ^3 X% rwindow1=hamming(100); %海明窗
window2=blackman(100); %blackman窗
4 k2 K) `8 J# M$ hnoverlap=20; %数据无重叠
* l& [6 x4 g6 `! Nrange='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);
8 Q/ J5 C* `0 g; {) Wplot_Pxx1=10*log10(Pxx1);
plot_Pxx2=10*log10(Pxx2);
9 Z% w. g3 C4 P5 B
$ R) t! Y# P, _7 l8 j8 Ffigure(1)
: ?  `$ j1 d4 J2 Y8 }  R) U2 bplot(f,plot_Pxx);
; L" b5 W7 j  p. U8 Y9 M% f& E, \
2 H: u- w1 Y6 `5 Xpause;
4 H0 Z" D5 v. |5 R8 a6 O
figure(2)
plot(f,plot_Pxx1);

pause;
7 n' J! [; Z, R+ w7 Y" _
figure(3)
plot(f,plot_Pxx2);

$ c- ]& W  O2 |& K  q

NingW

matlab 功率谱分析

Heaven_1

plot函数比较好