功率谱密度估计方法的matlab实现.doc
1 / 13 功率谱密度估计方法的MATLAB 实现 在应用数学和物理学中，谱密度、功率谱密度和能量谱密度是一个用于信号的通用概 念，它表示每赫兹的功率、每赫兹的能量这样的物理量纲。在物理学中，信号通常是波的 形式，例如电磁波、随机振动或者声波。当波的频谱密度乘以一个适当的系数后将得到每 单位频率波携带的功率，这被称为信号的功率谱密度(power spectral density, PSD)或者谱 功率分布(spectral power distribution, SPD) 。功率谱密度的单位通常用每赫兹的瓦特数 (W/Hz)表示，或者使用波长而不是频率，即每纳米的瓦特数(W/nm)来表示。信号的 功率谱密度当且仅当信号是广义的平稳过程的时候才存在。如果信号不是平稳过程，那么 自相关函数一定是两个变量的函数，这样就不存在功率谱密度，但是可以使用类似的技术 估计时变谱密度。信号功率谱的概念和应用是电子工程的基础，尤其是在电子通信系统中， 例如无线电和微波通信、雷达以及相关系统。因此学习如何进行功率谱密度估计十分重要， 借助于 Matlab 工具可以实现各种谱估计方法的模拟仿真并输出结果。下面对周期图法、修 正周期图法、最大熵法、Levinson 递推法和 Burg 法的功率谱密度估计方法进行程序设计及 仿真并给出仿真结果。 以下程序运行平台：Matlab Ra(8.5.0.197613) 一、周期图法谱估计程序 1、源程序 Fs=100000; %采样频率 100kHz N=1024; %数据长度 N=1024 n=0:N-1; t=n/Fs; xn=sin(2000*2*pi*t); %正弦波,f=2000Hz Y=awgn(xn,10); %加入信噪比为 10db的高斯白噪声 subplot(2,1,1); plot(n,Y) title( 信号 ) 2 / 13 xlabel( 时间 );ylabel( 幅度 ); grid on; window=boxcar(length(xn)); %矩形窗 nfft=N/4; %采样点数 [Pxx f]=periodogram(Y,window,nfft,Fs); %直接法 subplot(2,1,2); plot(f,10*log10(Pxx)); grid on; title([ 周期图法谱估计, ,int2str(N), 点 ]); xlabel( 频率(Hz) );ylabel( 功率谱密度 ); 2、仿真结果 时间 0 200 400 600 800 1000 1200 幅度 -2 -1 0 1 2 信号 频率(Hz) 10 4 0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5 功率谱密度 -80 -60 -40 -20 周期图法谱估计,1024点 二、修正周期图法(加窗)谱估计程序 1、源程序 Fs=100000; %采样频率 100kHz N=512; %数据长度 M=32; %汉明窗宽度 3 / 13 n=0:N-1; t=n/Fs; xn=sin(2000*2*pi*t); %正弦波,f=2000Hz Y=awgn(xn,10); %加入信噪比为 10db的高斯白噪声 subplot(2,1,1);subplot(2,1,1); plot(n,Y) title( 信号 ) xlabel( 时间 );ylabel( 幅度 ); grid on; window=hamming(M); %汉明窗 [Pxx f]=pwelch(Y,window,10,256,Fs); subplot(2,1,2); plot(f,10*log10(Pxx)); grid on; title([ 修正周期图法谱估计 N= ,int2str(N), M= ,int2str(M)]); xlabel( 频率(Hz) );ylabel( 功率谱密度 ); 2、仿真结果 4 / 13 时间 0 100 200 300 400 500 600 幅度 -2 -1 0 1 2 信号 频率(Hz) 10 4 0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5 功率谱密度 -60 -50 -40 -30 修正周期图法谱估计N=512 M=32 三、最大熵法谱估计程序 1、源程序 fs=1; %设采样频率 N=128; %数据长度 改变数据长度会导致分辨率的变化； f1=0.2*fs; %第一个sin信号的频率,f1/fs=0.2 f2=0.3*fs; %第二个sin信号的频率，f2/fs=0.2或者0.3 P=10; %滤波器阶数 n=1:N; s=sin(2*pi*f1*n/fs)+sin(2*pi*f2*n/fs); %s为原始信号 x=awgn(s,10); %x为观测信号，即对原始信号加入白噪声，信噪比10dB figure(1); %画出原始信号和观测信号 subplot(2,1,1); plot(s, b ),xlabel( 时间 ),ylabel( 幅度 ),title( 原始信号s ); grid; subplot(2,1,2); plot(x, r ),xlabel( 时间 ),ylabel( 幅度 ),title( 观测信号x ); 5 / 13 [Pxx1,f]=pmem(x,P,N,fs); %最大熵谱估计 figure(2); plot(f,10*log10(Pxx1)); xlabel( 频率(Hz) );ylabel( 功率谱(dB) ); title([ 最大熵法谱估计 模型阶数P= ,int2str(P), 数据长度 N= ,int2str(N)]); 2、仿真结果 时间 0 20 40 60 80 100 120 140 幅度 -2 -1 0 1 2 原始信号s 时间 0 20 40 60 80 100 120 140 幅度 -4 -2 0 2 4 观测信号x 6 / 13 频率(Hz) 0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 0.5 功率谱(dB) -15 -10 -5 0 5 10 15 最大熵法谱估计 模型阶数P=10 数据长度N=128 四、Levinson 递推法谱估计程序 1、源程序 fs=1; %设采样频率为 1 N=1000; %数据长度 改变数据长度会导致分辨率的变化； f1=0.2*fs; %第一个 sin 信号的频率,f1/fs=0.2 f2=0.3*fs; %第二个 sin 信号的频率，f1/fs=0.2 或者 0.3 M=16; %滤波器阶数的最大取值,超过则认为代价太大而放弃 L=2*N; %有限长
