matlab 功率谱分析

mutougeda

matlab 功率谱分析

1、直接法：
: U8 g4 @. p% V- e# {+ f; ]直接法又称周期图法，它是把随机序列x(n)的N个观测数据视为一能量有限的序列，直接计算x(n)的离散傅立叶变换，得X(k)，然后再取其幅值的平方，并除以N，作为序列x(n)真实功率谱的估计。
Matlab代码示例：
7 K7 Y& b/ _7 Y" l4 ^clear;
' a8 P( N( D' ?; j8 vFs=1000; %采样频率
% S9 @' B0 G) }- \n=0:1/Fs:1;
%产生含有噪声的序列
: Q' L6 \: P' o% ^0 ?5 j3 Vxn=cos(2*pi*40*n)+3*cos(2*pi*100*n)+randn(size(n));
0 H3 @. W9 R  Y$ c$ h# Swindow=boxcar(length(xn)); %矩形窗
nfft=1024;
[Pxx,f]=periodogram(xn,window,nfft,Fs); %直接法
plot(f,10*log10(Pxx));

2、间接法：
: Q$ q  \# V* t* j! d; q间接法先由序列x(n)估计出自相关函数R(n)，然后对R(n)进行傅立叶变换，便得到x(n)的功率谱估计。
5 V- B4 W; W; F% Z8 jMatlab代码示例：
9 L! r! C3 Z8 s# j! gclear;
Fs=1000; %采样频率
n=0:1/Fs:1;
' Z  o2 I! ~3 L3 G) u7 f: F%产生含有噪声的序列
8 r- N$ W/ W7 p/ ?2 i9 r) Zxn=cos(2*pi*40*n)+3*cos(2*pi*100*n)+randn(size(n));
2 o" \5 U5 i$ b- \3 }5 ?3 Wnfft=1024;
cxn=xcorr(xn,'unbiased'); %计算序列的自相关函数
7 z0 h4 E. Z# ?; c# U/ PCXk=fft(cxn,nfft);
+ y% @' N1 p! e! Q+ d; _0 N, qPxx=abs(CXk);
index=0:round(nfft/2-1);
k=index*Fs/nfft;
plot_Pxx=10*log10(Pxx(index+1));
plot(k,plot_Pxx);
6 Z% ^- q9 b) c- ^7 i
9 c6 K4 c3 X4 J, X" |3、改进的直接法：
. {2 B' n- D. O; m对于直接法的功率谱估计，当数据长度N太大时，谱曲线起伏加剧，若N太小，谱的分辨率又不好，因此需要改进。
) ?2 y8 O. S8 q7 y6 M5 Y3.1、Bartlett法
Bartlett平均周期图的方法是将N点的有限长序列x(n)分段求周期图再平均。
Matlab代码示例：
clear；
, p3 ^7 |  L* u' H9 q. {; GFs=1000;
n=0:1/Fs:1;
xn=cos(2*pi*40*n)+3*cos(2*pi*100*n)+randn(size(n));
; x8 R: P1 `8 i0 xnfft=1024;
window=boxcar(length(n)); %矩形窗
; d& ?) z) L0 U$ `! `: Xnoverlap=0; %数据无重叠
p=0.9; %置信概率
. I# M0 G, w3 w8 p/ U$ q[Pxx,Pxxc]=psd(xn,nfft,Fs,window,noverlap,p);
0 ~8 y" O' s5 h: I# i1 Oindex=0:round(nfft/2-1);
k=index*Fs/nfft;
, h& O  x0 k+ u+ iplot_Pxx=10*log10(Pxx(index+1));
5 V- ^6 F7 b7 V( B! X, V2 g. |plot_Pxxc=10*log10(Pxxc(index+1));
) N. b! w( {! }3 O/ cfigure(1)
+ F6 a) d% O- Gplot(k,plot_Pxx);
  b8 B% b$ e( @8 }) E  m0 `pause;
' ^1 g7 e0 y4 S+ B7 p- G) `/ Jfigure(2)
' y' j: B- q- D8 Vplot(k,[plot_Pxx plot_Pxx-plot_Pxxc plot_Pxx+plot_Pxxc]);

7 Z  J, v( m1 j" ~9 g2 i3.2、Welch法
Welch法对Bartlett法进行了两方面的修正，一是选择适当的窗函数w(n)，并再周期图计算前直接加进去，加窗的优点是无论什么样的窗函数均可使谱估计非负。二是在分段时，可使各段之间有重叠，这样会使方差减小。
Matlab代码示例：
% O0 D6 N# m: @6 L( N$ Xclear;
Fs=1000;
$ o4 b( v! z9 zn=0:1/Fs:1;
xn=cos(2*pi*40*n)+3*cos(2*pi*100*n)+randn(size(n));
* g+ o+ F9 G( @& ?0 S- l3 b! infft=1024;
window=boxcar(100); %矩形窗
( W0 W8 X# `; ^window1=hamming(100); %海明窗
window2=blackman(100); %blackman窗
noverlap=20; %数据无重叠
/ O/ G$ o$ a7 c& S/ [  Y- k; [6 Drange='half'; %频率间隔为[0 Fs/2]，只计算一半的频率
5 A! T9 ]: ]9 V4 B4 t7 u0 T[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);
$ x! r# w$ b4 m& N3 Nplot_Pxx=10*log10(Pxx);
plot_Pxx1=10*log10(Pxx1);
plot_Pxx2=10*log10(Pxx2);

5 O7 V4 H. i8 j, u8 Efigure(1)
7 [& v3 T% I# p+ n/ ]: vplot(f,plot_Pxx);
) ~7 B4 T( r! d% c
pause;
+ V% r; f- [7 ~, r; `- }- r; {" m
figure(2)
plot(f,plot_Pxx1);
5 W* |, I, [! S
2 O- @6 t) t1 [7 S( v; d7 _# Cpause;

figure(3)
% @% T# \! X1 Fplot(f,plot_Pxx2);

NingW

matlab 功率谱分析

Heaven_1

plot函数比较好