matlab 功率谱分析

mutougeda

matlab 功率谱分析
& q/ `3 Z: C0 C

1、直接法：
直接法又称周期图法，它是把随机序列x(n)的N个观测数据视为一能量有限的序列，直接计算x(n)的离散傅立叶变换，得X(k)，然后再取其幅值的平方，并除以N，作为序列x(n)真实功率谱的估计。
& z, u% h+ J: G( c' a% FMatlab代码示例：
8 M' Q# D( X2 `9 X1 ^3 o9 @6 N. x' aclear;
, [& N+ A# c" N# s, cFs=1000; %采样频率
, [/ Y" L# D$ o% `/ f: g. |! V3 wn=0:1/Fs:1;
%产生含有噪声的序列
xn=cos(2*pi*40*n)+3*cos(2*pi*100*n)+randn(size(n));
0 l! y5 c9 R, Z- @, pwindow=boxcar(length(xn)); %矩形窗
nfft=1024;
& N) W. }* L0 a; {[Pxx,f]=periodogram(xn,window,nfft,Fs); %直接法
% v. x1 x/ X) L2 D; s* dplot(f,10*log10(Pxx));

2 @( U1 P2 X  k2 R- M0 }2、间接法：
' l/ H1 h- L; m/ ]/ P间接法先由序列x(n)估计出自相关函数R(n)，然后对R(n)进行傅立叶变换，便得到x(n)的功率谱估计。
Matlab代码示例：
' ]. w/ j: Y3 O1 |. `: n# eclear;
Fs=1000; %采样频率
n=0:1/Fs:1;
; n4 D5 k5 U& b" g0 b%产生含有噪声的序列
0 P" D1 X1 t+ X) r$ ]- exn=cos(2*pi*40*n)+3*cos(2*pi*100*n)+randn(size(n));
nfft=1024;
6 P; J* U9 O1 Z4 Z1 ?cxn=xcorr(xn,'unbiased'); %计算序列的自相关函数
0 b* c" c' h5 @CXk=fft(cxn,nfft);
* m, e; D/ \$ x& W% qPxx=abs(CXk);
. g- q; J/ I/ K& iindex=0:round(nfft/2-1);
k=index*Fs/nfft;
3 X& a6 {  d+ ~# i% a: F9 u* Gplot_Pxx=10*log10(Pxx(index+1));
# C5 y/ u1 N7 u. j! h) Z9 E6 zplot(k,plot_Pxx);
2 g( V0 {2 R/ x! j& i
* S. g- Z* R: Z' ^3、改进的直接法：
4 w: A! h8 y: j对于直接法的功率谱估计，当数据长度N太大时，谱曲线起伏加剧，若N太小，谱的分辨率又不好，因此需要改进。
3.1、Bartlett法
* o- Y5 R+ M3 V3 a3 C" Y' HBartlett平均周期图的方法是将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));
" U& j: ^  u8 R* b: }  Onfft=1024;
window=boxcar(length(n)); %矩形窗
& A3 k& y2 E* c; t/ v$ X: wnoverlap=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));
& H) J; t  x# {( v0 Wfigure(1)
plot(k,plot_Pxx);
pause;
6 e# t. V7 m  lfigure(2)
- Z7 F% e' l$ o1 ~2 L( [plot(k,[plot_Pxx plot_Pxx-plot_Pxxc plot_Pxx+plot_Pxxc]);
+ t) W4 V" F/ m. x5 j
3.2、Welch法
/ S' A6 m6 t! ]# E( q8 NWelch法对Bartlett法进行了两方面的修正，一是选择适当的窗函数w(n)，并再周期图计算前直接加进去，加窗的优点是无论什么样的窗函数均可使谱估计非负。二是在分段时，可使各段之间有重叠，这样会使方差减小。
& t) \4 ~* [/ b) o+ cMatlab代码示例：
clear;
$ m: x, h: m- U9 B6 z6 SFs=1000;
n=0:1/Fs:1;
' e' ]; J% ]! _3 f+ @xn=cos(2*pi*40*n)+3*cos(2*pi*100*n)+randn(size(n));
4 D' f, _$ D9 O. |$ _4 anfft=1024;
window=boxcar(100); %矩形窗
window1=hamming(100); %海明窗
$ g6 Q. `& N5 }# N+ cwindow2=blackman(100); %blackman窗
noverlap=20; %数据无重叠
% t6 z$ r( G* v3 m1 z. d/ yrange='half'; %频率间隔为[0 Fs/2]，只计算一半的频率
, }) l3 i, c2 }[Pxx,f]=pwelch(xn,window,noverlap,nfft,Fs,range);
[Pxx1,f]=pwelch(xn,window1,noverlap,nfft,Fs,range);
; H- `* q( {  G7 [  I  q[Pxx2,f]=pwelch(xn,window2,noverlap,nfft,Fs,range);
, u) z1 c% n  z$ ~1 Kplot_Pxx=10*log10(Pxx);
plot_Pxx1=10*log10(Pxx1);
) S1 a' i' ]! u0 O  pplot_Pxx2=10*log10(Pxx2);

3 C4 A, H$ z  Hfigure(1)
plot(f,plot_Pxx);

pause;
& b& L6 [2 \6 X( u8 O- u4 K
) q1 _1 s, G+ r, m) Vfigure(2)
# G/ S9 }. n0 Q3 ~plot(f,plot_Pxx1);

" H5 H7 A! f% Z6 L1 Z2 c" \' }pause;
+ H; @  _. m% p" e! K# O# i
figure(3)
plot(f,plot_Pxx2);

7 }+ A# V1 ?+ u1 x
+ J! Y6 V) {4 f+ [! Y+ k# A

NingW

matlab 功率谱分析

Heaven_1

plot函数比较好