Burg算法
基于BURG算法的谱估计研究及其MATLAB实现课件
三、经典谱估计
当f1=20, f2=40时得到的仿真结果如图3.5所示。 当f1=100, f2=1000时得到的仿真结果如图3.6 所示。
图3.5
图3.6
三、经典谱估计
当噪声方差 =1和 =10的仿真结果分别如图 3.7所示和图3.8所示。
图3.7
图3.8
三、经典谱估计
仿真结果: 间接法实现时在噪声信号很小的情况 下, 图3.5得到的谱线基本能分辨出两个频率值来, 但是也出现大量假峰。图3.6也能分辨出两个频率值, 假峰减少, 说明间接法在这种情况下得到的效果要
声信号增大到原来的100倍时就无法分辨两个频率 值, 而且通过多次仿真看出, 当噪声信号增大到10倍 以
后就不能分辨出两个频率点的峰值。
三、经典谱估计
■ 3.2 间接法及MATLAB仿真
结果
的估
■ 间接法又称自相关法, 记 为对计, 即
当M较小时, 上式的计算量不是很大, 因此, 此 方法是在FFT问世之前（即周期图被广泛应用之前） 常用的谱估计方法。
四、现代谱估计
当阶数=10时得到的仿真结果如图4.1所示。 当阶数=15时得到的仿真结果如图4.2所示。
图4.1
图4.2
四、现代谱估计
仿真结果: Levinson-Durbin算法得到的谱线波 动性小, 能很好的分辨出两个频率值, 而且没有出 现假峰现象。当增大阶数时得到的结果跟阶数小的 结果不相上下, 并没有增大频率分辨率, 反而增大 了计算次数。
仿真结果比较可得: 四种算法都能分辨出两个频 率值，可以明显看出直接法和间接法得到的仿真结果 出现了大量假峰，且谱的波动性较大，而 LevinsonDurbin算法和BURG算法得到的仿真结果很平滑，没有 出现假峰现象，能清楚的分辨出两个频率点的值，分 辨率远比经典谱估计要好，所以从仿真的结果也可以 看出现代谱估计性能比经典谱估计好。
现代信号处理大作业题目+答案
研究生“现代信号处理”课程大型作业（以下四个题目任选三题做）1. 请用多层感知器（MLP ）神经网络误差反向传播（BP ）算法实现异或问题（输入为[00;01;10;11]X T =，要求可以判别输出为0或1），并画出学习曲线。
其中，非线性函数采用S 型Logistic 函数。
2. 试用奇阶互补法设计两带滤波器组(高、低通互补)，进而实现四带滤波器组；并画出其频响。
滤波器设计参数为：F p ＝1.7KHz , F r ＝2.3KHz , F s ＝8KHz , A rmin ≥70dB 。
3. 根据《现代数字信号处理》（姚天任等，华中理工大学出版社，2001）第四章附录提供的数据(pp.352-353)，试用如下方法估计其功率谱，并画出不同参数情况下的功率谱曲线：1） Levinson 算法2） Burg 算法3） ARMA 模型法4） MUSIC 算法4. 图1为均衡带限信号所引起失真的横向或格型自适应均衡器(其中横向FIR 系统长M =11), 系统输入是取值为±1的随机序列)(n x ，其均值为零；参考信号)7()(-=n x n d ；信道具有脉冲响应：12(2)[1cos()]1,2,3()20 n n h n W π-⎧+=⎪=⎨⎪⎩其它式中W 用来控制信道的幅度失真(W = 2～4, 如取W = 2.9,3.1,3.3,3.5等)，且信道受到均值为零、方差001.02=v σ(相当于信噪比为30dB)的高斯白噪声)(n v 的干扰。
试比较基于下列几种算法的自适应均衡器在不同信道失真、不同噪声干扰下的收敛情况(对应于每一种情况,在同一坐标下画出其学习曲线):1） 横向/格-梯型结构LMS 算法2） 横向/格-梯型结构RLS 算法并分析其结果。
图1 横向或格-梯型自适应均衡器参考文献[1] 姚天任, 孙洪. 现代数字信号处理[M]. 武汉: 华中理工大学出版社, 2001[2] 杨绿溪. 现代数字信号处理[M]. 北京: 科学出版社, 2007[3] S. K. Mitra. 孙洪等译. 数字信号处理——基于计算机的方法(第三版)[M]. 北京: 电子工业出版社, 2006[4] S.Haykin, 郑宝玉等译. 自适应滤波器原理(第四版)[M].北京: 电子工业出版社, 2003[5] J. G. Proakis, C. M. Rader, F. Y. Ling, etc. Algorithms for Statistical Signal Processing [M].Beijing: Tsinghua University Press, 2003一、请用多层感知器（MLP）神经网络误差反向传播（BP）算法实现异或问题（输入为[00;01;10;11]，要求可以判别输出为0或1），并画出学习曲线。
基于Burg算法的最大熵谱估计
基于Burg 算法的最大熵谱估计一、 实验目的使用Matlab 平台实现基于Burg 算法的最大熵谱估计二、 Burg 算法原理现代谱估计是针对经典谱估计方差性能较差、分辨率较低的缺点提出并逐渐发展起来的，其分为参数模型谱估计和非参数模型谱估计。
而参数模型谱估计主要有AR 模型、MA 模型、ARMA 模型等，其中AR 模型应用最多。
ARMA 模型功率谱的数学表达式为：212121/1)(∑∑=-=-++=p i i j i q i i j i j e a e b e P ωωωσ其中，P(e j ω)为功率谱密度；s 2是激励白噪声的方差；a i 和b i 为模型参数。
若ARMA 模型中b i 全为0，就变成了AR 模型，又称线性自回归模型，其是一个全极点模型： 2121/)(∑=-+=p i i j i j e a e P ωωσ研究表明，ARMA 模型和MA 模型均可用无限阶的AR 模型来表示。
且AR 模型的参数估计计算相对简单。
同时，实际的物理系统通常是全极点系统。
要利用AR 模型进行功率谱估计，必须由Yule - Walker 方程求得AR 模型的参数。
而目前求解Yule - Walker 方程主要有三种方法: Levinson-Durbin 递推算法、Burg 算法和协方差方法。
其中Burg 算法计算结果较为准确，且对于短的时间序列仍能得到较正确的估计，因此应用广泛。
研究最大熵谱估计时，Levinson 递推一直受制于反射系数K m 的求出。
而Burg 算法秉着使前、后向预测误差平均功率最小的基本思想，不直接估计AR 模型的参数，而是先估计反射系数K m ，再利用Levinson 关系式求得AR 模型的参数，继而得到功率谱估计。
Burg 定义m 阶前、后向预测误差为：∑=-=mi m m i n x i a n f 0)()()( (1)∑=*--=mi m i n x i m a n g m 0)()()( (2) 由式（1）和（2）又可得到前、后预测误差的阶数递推公式：)1()()(11-+=--n g K n f n f m m m m (3))1()()(11-+=--*n g n f K n g m m m m (4)定义m 阶前、后向预测误差平均功率为：∑=+=Nmn m m m n g n f P ])()([2122(5) 将阶数递推公式（3）和（4）代入（5），并令0=∂∂mmK P ，可得∑∑+=--+=*---+--=N m n m m Nm n m m m n g n f n g n f K 12121111])1()([21)1()((6)三、 Burg 算法递推步骤Burg 算法的具体实现步骤：步骤1 计算预测误差功率的初始值和前、后向预测误差的初始值，并令m = 1。
Y-W方程的L-D递推算法和伯格递推算法的比较
(1) 确定初始条件
e [ k ] e [ k ] y[ k ]
f 0 b 0
2 0
1 N
N 1 k 0
y 2 [k ]
(2) 从p=1开始迭代计算: 计算AR模型参数
Kp 2 e fp 1 [ k ]e b p 1 [ k 1]
谱估计结果——p=40,N=128
Periodogram 60 40 20 0 -20 -40 -60 0 0.1 0.2 0.3 0.4 0.5 Yule 60 40 20 0 -20 -40 -60 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 0.6 0.7 0.8 0.9 1
[Pxx,f] = pburg (x,p,NFFT,Fs) x：进行功率谱估计的输入有限长序列； p: 模型的阶数 NFFT：DFT的点数； Fs ：绘制功率谱曲线的抽样频率，默认值为1； Pxx：功率谱估计值； f：Pxx值所对应的频率点
利用Burg法进行谱估计程序
N=512;NFFT=1024;Fs=2;p=40; n=0:N-1;randn('state',0); x=cos(0.3*pi*n)+cos(0.32*pi*n)+randn(size(n)); [P,f]=pyulear(x,p,NFFT,2); [Pw, f2]=pburg(x,p,NFFT,2); subplot(211);plot(f,10*log(P));grid;title('L-D'); axis([0 1 -30 60]); subplot(212);plot(f2,10*log(Pw));grid;title('Burg '); axis([0 1 -30 60]);
Burg算法
功率谱估计的古典算法与现代算法的比较——选取周期图法与Burg算法为例现代信号分析中, 对于常见的具有各态历经的平稳随机信号, 不可能用清楚的数学关系式来描述, 但可以利用给定的 N 个样本数据估计一个平稳随机信号的功率谱密度叫做功率谱估计(PSD)。
功率谱估计可以分为经典功率谱估计(非参数估计)和现代功率谱估计(参数估计)。
一、古典功率谱估计古典功率谱估计是将数据工作区外的未知数据假设为零, 相当于数据加窗经典功率谱估计方法分为: 相关函数法(BT 法)、周期图法以及两种改进的周期图估计法。
1、相关法相关法是以相关函数为媒介来计算功率谱的，所以又叫间接法，它的理论基础是维纳--辛钦定理。
先对数据工作区外的未知数据赋值为零，再由序列x(n)估计出自相关函数R(n)，最后对R(n)进行傅立叶变换, 便得到 x(n)的功率谱估计。
2、周期图法周期图法是由获得的N点数据构成的有限长序列直接求fft得其频谱，取频谱幅度的平方再除以N，以此作为对x(n)真实功率谱的估计。
3、改进的周期图法改进的周期图法的主要途径是平滑和平均。
平滑是用一个适当的窗函数与算出的功率谱进行卷积，使谱线平滑，这种方法得出的谱估计是无偏的，方差也小，但分辨率下降；平均就是将截取的数据段再分成L个平均的小段，分别计算功率谱后取功率谱的平均，当L趋于无穷大的时候，L个平均的方差趋于零，可以达到一致谱估计的目的。
由于存在旁瓣，会产生两个后果：一是功率谱主瓣能量泄露到旁瓣使谱估计的方差增大，二是与旁瓣卷积后得到的功率谱完全属于干扰，严重情况下，强信号与旁瓣的卷积可能大于弱信号与主瓣的卷积，使弱信号淹没在强信号的干扰中无法检测出来。
这是古典法谱估计的主要缺点，即便是改进的周期图法也无法克服分辨率低的缺点。
我们从中选取周期图法作比较，其算法实现如下：Fs=600; %采样频率n=0:1/Fs:1;%产生含有噪声的序列xn=cos(2*pi*40*n)+cos(2*pi*90*n)+0.1*randn(size(n));n=1:length(xn);figure(1);subplot(2,1,1);plot(n,xn);window=boxcar(length(xn));%矩形窗nfft=1024;[Pxx,f]=periodogram(xn,window,nfft,Fs);subplot(2,1,2);plot(f,10*log10(Pxx));得到的图形为：二、现代谱估计参数模型法是现代谱估计中的主要内容，AR 模型参数的求解有三种方法：自相关法、Burg 递推算法和改进协方差法。
现代信号处理大作业
现代信号处理大型作业一．试用奇阶互补法设计两带滤波器组(高、低通互补)，进而实现四带滤波器组；并画出其频响。
滤波器设计参数为：F p ＝1.7KHz , F r ＝2.3KHz , F s ＝8KHz , A rmin ≥70dB 。
（一）、分析与通常的滤波器相比，互补滤波器具有优良的结构特性和结构特性，具有较低的噪声能量和系数敏感性，其定义如下：一组滤波器H 12(),(),.......()Z H Z H Z n 如果满足下式：He Kjw k n(),==∑110<w<2π 则称这组滤波器为幅度互补滤波器;如果满足下式：He kjw k n()=∑=121, 0<w<2π则称这组滤波器为功率互补滤波器，同时互补滤波器还应该满足：Hz A z kk n()()=∑=1其中A(z)为全通函数，适当的选择全通函数，可以使两带函数具有所需要的低通和高通特性。
（二）、设计步骤(1) 对Fp 、Fr 进行预畸);();(''FsFrtg FsFptg r p ∏=Ω∏=Ω(2) 计算'''*r p c ΩΩ=Ω，判断'c Ω是否等于1，即该互补滤波器是否为互补镜像滤波器(3) 计算相关系数⎪⎩⎪⎨⎧-==+++=+-=-=ΩΩ=--=偶数)N 为(;21奇数)N 为 (;;lg /)16/1lg(;150152;1121;1;;])110)(110[(1213090500''02'''211－min1.0min1.0i i u q k N q q q q q k k q k k k k rp Ar Ap;)2cos()1(21))12(sin()1(21)1(21'2∑∑∞=∞=+-++-=Ωm mm m m m m i u Nm q u Nm q q ππ;42⎥⎦⎤⎢⎣⎡=N N;221N N N -⎥⎦⎤⎢⎣⎡=;)/1)(1(2'2'k k v i i i Ω-Ω-=12'1212,1;12N i v i i i =Ω+=--α 22'22,1;12N i v iii =Ω+=β (4) 互补镜像滤波器的数字实现;22i ii A αα+-=;22iii B ββ+-=1221,1;1)(N i ZA Z A Z H i i i =++=∏--22212,1;1)(N i ZB Z B Z Z H i i i =++=∏--- )];()([21)(21Z H Z H Z H L +=（三）、程序与结果 1． 二带滤波器组 （1） 源程序： clear; clf;Fp=1700;Fr=2300;Fs=8000; Wp=tan(pi*Fp/Fs); Wr=tan(pi*Fr/Fs); Wc=sqrt(Wp*Wr); k=Wp/Wr;k1=sqrt(sqrt(1-k^2)); q0=0.5*(1-k1)/(1+k1);q=q0+2*q0^5+15*q0^9+150*q0^13; N=11;N2=fix(N/4); M=fix(N/2); N1=M-N2; for jj=1:M a=0;for m=0:5a=a+(-1)^m*q^(m*(m+1))*sin((2*m+1)*pi*jj/N);%N is odd, u=j end ab=0;for m=1:5b=b+(-1)^m*q^(m^2)*cos(2*m*pi*jj/N); end bW(jj)=2*q^0.25*a/(1+2*b);V(jj)=sqrt((1-k*W(jj)^2)*(1-W(jj)^2/k)); endfor i=1:N1alpha(i)=2*V(2*i-1)/(1+W(2*i-1)^2); endfor i=1:N2beta(i)=2*V(2*i)/(1+W(2*i)^2); endfor i=1:N1a(i)=(1-alpha(i)*Wc+Wc^2)/(1+alpha(i)*Wc+Wc^2); endfor i=1:N2b(i)=(1-beta(i)*Wc+Wc^2)/(1+beta(i)*Wc+Wc^2); endw=0:0.0001:0.5;LP=zeros(size(w));HP=zeros(size(w));for n=1:length(w)z=exp(j*w(n)*2*pi);H1=1;for i=1:N1H1=H1*(a(i)+z^(-2))/(1+a(i)*z^(-2)) ;endH2=1/z;for i=1:N2H2=H2*(b(i)+z^(-2))/(1+b(i)*z^(-2));endLP(n)=abs((H1+H2)/2);HP(n)=abs((H1-H2)/2);endplot(w,LP,'b',w,HP,'r');hold on;xlabel('digital frequency');ylabel('amptitude');（2）运行结果：见图1图1 二带数字滤波器组2．四带滤波器组（1）源程序：clf;Fp=1700;Fr=2300;Fs=8000;Wp=tan(pi*Fp/Fs);Wr=tan(pi*Fr/Fs);Wc=sqrt(Wp*Wr);k=Wp/Wr;k1=sqrt(sqrt(1-k^2));q0=0.5*(1-k1)/(1+k1);q=q0+2*q0^5+15*q0^9+150*q0^13;N=11;N2=fix(N/4);M=fix(N/2);N1=M-N2;for jj=1:Ma=0;for m=0:5a=a+(-1)^m*q^(m*(m+1))*sin((2*m+1)*pi*jj/N); % N is odd, u=jendb=0;for m=1:5b=b+(-1)^m*q^(m^2)*cos(2*m*pi*jj/N);endW(jj)=2*q^0.25*a/(1+2*b);V(jj)=sqrt((1-k*W(jj)^2)*(1-W(jj)^2/k));Endfor i=1:N1alpha(i)=2*V(2*i-1)/(1+W(2*i-1)^2);endfor i=1:N2beta(i)=2*V(2*i)/(1+W(2*i)^2);endfor i=1:N1a(i)=(1-alpha(i)*Wc+Wc^2)/(1+alpha(i)*Wc+Wc^2);endfor i=1:N2b(i)=(1-beta(i)*Wc+Wc^2)/(1+beta(i)*Wc+Wc^2);endw=0:0.0001:0.5;LLP=zeros(size(w));LHP=zeros(size(w));HLP=zeros(size(w));HHP=zeros(size(w));for n=1:length(w)z=exp(j*w(n)*2*pi);H1=1;for i=1:N1H1=H1*(a(i)+z^(-2))/(1+a(i)*z^(-2)) ;endH21=1;for i=1:N1H21=H21*(a(i)+z^(-4))/(1+a(i)*z^(-4)) ;H2=1/z;for i=1:N2H2=H2*(b(i)+z^(-2))/(1+b(i)*z^(-2));endH22=1/(z^2);for i=1:N2H22=H22*(b(i)+z^(-4))/(1+b(i)*z^(-4));endLP=((H1+H2)/2);HP=((H1-H2)/2);LLP(n)=abs((H21+H22)/2*LP);LHP(n)=abs((H21-H22)/2*LP);HHP(n)=abs((H21+H22)/2*HP);HLP(n)=abs((H21-H22)/2*HP);endplot(w,LLP,'b',w,LHP,'r',w,HLP,'k',w,HHP,'m')hold onxlabel('digital frequency');ylabel('amptitude');（2）运行结果：见图2图2 四带数字滤波器组二、根据《现代数字信号处理》第四章提供的数据，试用如下方法估计其功率谱，并画出不同参数情况下的功率谱曲线：1）Levison算法2）Burg算法3） ARMA 模型法 4） MUSIC 算法 1 Levinson 算法Levinson 算法用于求解Yule-Walker 方程，是一种按阶次进行递推的算法，即首先以AR （0）和AR （1）模型参数作为初始条件，计算AR （2）模型参数；然后根据这些参数计算AR （3）参数，等等，一直到计算出AR （p ）模型参数为止，需要的运算量数量级为2p ，其中p 为AR 模型的阶数。
Yule-Walker法和Burg法对信号进行谱估计的Matlab仿真
摘要：用现代谱估计中的AR 模型参数法中的Y ule-Walker 法和Burg 法对信号进行谱估计，并用MATLAB 进行信号的仿真，对于不同阶数下的信号进行对比分析。
关键词：谱估计；AR 模型参数法；Yule-Walker 法；Burg 法；MA TLAB （一）原理 1.Y ule-Walker 法：将一平稳随机信号x(n)表示成一个白噪声w(n)激励一个因果稳定的可逆系统H(z)产生的输出，再由已知的x(n)及其自相关函数()x R m 来估计H(z)的参数，由（1）式，可以用H （z ）的参数来表示x(n)的功率谱。
22()|()|jw jw x S e H e σ= （1） AR 模型又称为自回归模型，系统函数H （z ）只有极点没有零点,P 阶AR 模型的系统函数为：1()1pii i G H z a z -==+∑在白噪声激励下的输出：1()()()pi i x n a x n i Gw n ==--+∑令预测误差为e(n):1()()()()pi i e n Gw n x n a x n i ===+-∑尤勒-沃克方程：121(1),1,2,...,()(),0pi x i x pi x i a R m m p R m a R i G m ==⎧--=⎪⎪=⎨⎪-+=⎪⎩∑∑可表示为下面矩阵式形式：2121(0)(1)(2)()(1)(0)(1)(1)0(2)(1)(0)(2)0()(1)(2)(0)0x x x x x x x x x x x x p x x x x R R R R p a R R R R p a R R R R p a R p R p R p R σ⎡⎤⎡⎤⎡⎤⎢⎥⎢⎥⎢⎥-⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥-=⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥--⎣⎦⎣⎦⎣⎦只要已知或估计出p+1个自相关函数就可以由此方程解出p+1个模型参数{}21,2,,aa σ 。
Levinson-Durbin 递推算法是解尤勒-沃克方程的快速有效的算法，这种算法利用方程组系数矩阵所具有的一系列好的性质，使运算量大大减少。
(实验六 随机信号功率谱分析)
实验报告实验课程：数字信号处理实验开课时间：2020—2021 学年秋季学期实验名称：随机信号功率谱分析实验时间： 2020年9月30日星期三学院：物理与电子信息学院年级：大三班级：182 学号：1843202000234 姓名：武建璋一、实验预习实验目的要求深刻理解随机信号的特性，掌握随机信号功率谱估计的基本原理，灵活运用各种随机信号功率谱估计的基本方法。
实验仪器用具装有Matlab的计算机一台实验原理功率谱估计是随机信号处理中的一个重要的研究和应用领域.功率谱估计基本上可以非参数估计的经典方法和参数估计的近代方法.典型功率谱估计是基于FFT 算法的非参数估计,对足够长的记录数据效果较好。
在工程实际中，经典功率谱估计法获得广泛应用的是修正期图发。
该方法采取数据加窗处理再求平均的办法。
通过求各段功率谱平均，最后得到功率谱计P（m），即：式中：为窗口函数ω[k]的方差。
K表示有重叠的分数段。
由于采用分段加窗求功率谱平均，有效地减少了方差和偏差，提高了估计质量，使修正周期图法在经典法中得到普遍应用。
但在估计过程存在两个与实际不符的假设，即（1）利用有限的N个观察数据进行自相关估计，隐含着在已知N个数据之外的全部数据均为零的假设。
（2）假定数据是由N个观察数据以N为周期的周期性延拓。
同时在计算过程中采用加窗处理，使得估计的方差和功率泄露较大，频率分辨率较低，不适用于短系列的谱分析和对微弱信号的检测。
近代谱估计是建立在随机信号参数模型的基础上，通过信号参数模型或预测误差滤波器（一步预测器）参数的估计，实现功率谱估计。
由于既不需要加窗，又不需要对相关函数的估计进行如经典法那样的假设，从而减少公里泄露，提高了频谱分辨率。
常用的参数模型有自回归（AR）模型、滑动平均（MA）模型、自回归滑动平均（ARMA）模型。
其中AR模型是基本模型，求解AR模型的参数主要有L—D算法和Burg算法。
1.某随机信号由两余弦信号与噪声构成x(t)=cos(20*pi*t)+cos(40*pi*t)+s(t)式中：s(t)是均值为0、方差为1的高斯白噪声。
信号分析与处理ch1_6
power spectrum (dB)
power spectrum (dB)
实验六 随机信号功率谱分析
二、实验原理256 sample section,128 sample overlap 30 20
实验六 随机信号功率谱分析
二、实验原理
MATLAB也提供了近代功率谱估计的函数。
Burg算法可以利用函数arburg实现，其调用格式为
[A,E,K]= ARBURG(x, ORDER)
其中： ORDER为AR模型的阶数；
E为预测误差；K为反射系数。
x为观测序列；
实验六 随机信号功率谱分析
三、实验内容
Averwag=edhmaondniifinegd(p2e5ri6o)d';ogram,256 sample section,128sample overlap,hanning window Px3x03 = (abs(fft(w.*x( 1:256))).^2 + abs(fft(w.*x(129:384))).^2 + 20 abs(fft(w.*x(257:512))).^2 + abs(fft(w.*x(385:640))).^2 + 10 abs(fft(w.*x(513:768))).^2 + abs(fft(w.*x(641:896))).^2 ) ； Pxx03=Pxx3/ (norm(w)^2*6); su-1b0plot(2,1,2); pl-o20t([0:255]*Fs/256,10*log10(Pxx3)); axis([0,1023,-30,30]);
(完整版)功率谱估计性能分析及Matlab仿真
功率谱估计性能分析及Matlab 仿真1 引言随机信号在时域上是无限长的，在测量样本上也是无穷多的，因此随机信号的能量是无限的，应该用功率信号来描述。
然而，功率信号不满足傅里叶变换的狄里克雷绝对可积的条件，因此严格意义上随机信号的傅里叶变换是不存在的。
因此，要实现随机信号的频域分析，不能简单从频谱的概念出发进行研究，而是功率谱[1]。
信号的功率谱密度描述随机信号的功率在频域随频率的分布。
利用给定的N 个样本数据估计一个平稳随机信号的功率谱密度叫做谱估计。
谱估计方法分为两大类：经典谱估计和现代谱估计。
经典功率谱估计如周期图法、自相关法等，其主要缺陷是描述功率谱波动的数字特征方差性能较差，频率分辨率低。
方差性能差的原因是无法获得按功率谱密度定义中求均值和求极限的运算[2]。
分辨率低的原因是在周期图法中，假定延迟窗以外的自相关函数全为0。
这是不符合实际情况的，因而产生了较差的频率分辨率。
而现代谱估计的目标都是旨在改善谱估计的分辨率，如自相关法和Burg 法等。
2 经典功率谱估计经典功率谱估计是截取较长的数据链中的一段作为工作区，而工作区之外的数据假设为0，这样就相当将数据加一窗函数，根据截取的N 个样本数据估计出其功率谱[1]。
2.1 周期图法( Periodogram )Schuster 首先提出周期图法。
周期图法是根据各态历经的随机过程功率谱的定义进行的谱估计。
取平稳随机信号()x n 的有限个观察值(0),(1),...,(1)x x x n -，求出其傅里叶变换10()()N j j n N n X e x n e ωω---==∑然后进行谱估计21()()j N S X e Nωω-= 周期图法应用比较广泛，主要是由于它与序列的频谱有直接的对应关系，并且可以采用FFT 快速算法来计算。
但是，这种方法需要对无限长的平稳随机序列进行截断，相当于对其加矩形窗，使之成为有限长数据。
同时，这也意味着对自相关函数加三角窗，使功率谱与窗函数卷积，从而产生频谱泄露，容易使弱信号的主瓣被强信号的旁瓣所淹没，造成频谱的模糊和失真，使得谱分辨率较低[1]。
谱估计(现代)
ak xx (m k ) Ex(n) (n m)
k 1
p
而
m0 0, E x(n) (n m) 2 , m 0
•Yule-Walker方程的推导
故
p a k xx (m k ) , m 0 k 1 xx (m) p a (k ) 2 , m 0 k xx k 1 或
p
2
需要推导AR参数与 xx (m)之间的关系。
3.1
• 估计方法
自回归模型法
2 与xx (m)乊间的关系 参数a1, a2, a3, …, ap及 ——Yule-Walker方程
已知：自相关函数 已知： 自相关函数
Yule-Walker方程
要求： AR模型的阶数p，以及p个AR 要求： AR模型的阶数p，以及p个 AR 参数a(i),激励源方差 2 参数a(k),激励源方差
3.2
最大熵谱估计法
• 基本思想——熵
代表一种不定度； 最大熵为最大不定度，即它的时间序列最随机， 它的PSD应是最平伏(最白色)。 Shannon对熵的定义： 当x的取值为离散的时，熵H定义为
H pi ln pi
i
pi：出现状态i 的概率。
当x的取值为连续的时，熵H定义为
p(x)：概率密度 函数
(n)
...
z-1 a1
z-1
z-1
a2
...
ap
3.1
自回归模型法
q
• MA(Moving Average)模型 ——全零点模型
x(n) bl (n l )
l 0
H ( z ) B( z ) 1 bl z k
matlab实验之求均值-方差
实验一 随机信号的数字特征分析一、 实验目的1.了解随机信号自身的特性，包括均值（数学期望）、方差、均方值等；2. 掌握随机信号的分析方法；二、实验原理1.均值测量方法均值ˆx m表示集合平均值或数学期望值。
基于随机过程的各态历经性，最常用的方法是取N 个样本数据并简单地进行平均，即101ˆ[]N x d i m X i N-==∑ 其中，样本信号的采样数据记为[](,)d X i X iT ξ=，s T 为采样间隔。
2.均方误差的测量方法随机序列的均方误差定义为： 2211()lim ()N i N i E X x n N →∞==∑ 3.方差测量方法如果信号的均值是已知的，则其方差估计设计为12201ˆ([])N x X d i X i m N σ-==-∑ 它是无偏的与渐进一致的。
三、实验内容利用MATLAB 中的伪随机序列产生函数randn()产生多段1000点的序列，编制一个程序，计算随机信号的数字特征，包括均值、方差、均方值、最后把计算结果平均，绘制数字特征图形。
源程序如下：clear all;clc;%产生50个1000以内点的伪随机序列x=randn(50,1000);%计算随机产生的50个点序列的均值，方差，均方average=zeros(1,50);variance=zeros(1,50);square=zeros(1,50);%计算均值for i=1:50for j=1:1000average(i)=average(i)+x(i,j);endaverage(i)=average(i)/1000;end%计算方差for i=1:50for j=1:1000variance(i)=variance(i)+(x(i,j)-average(i)).^2; endvariance(i)=variance(i)/1000;end%计算均方值for i=1:50for j=1:1000square(i)=square(i)+x(i,j).^2;endsquare(i)=square(i)/1000;endEX=sum(average)/50;DX=sum(variance)/50;RMS=sum(square)/50;plot(average);title('50个随机序列的均值');figure;plot(variance);title('50个随机序列的方差');figure;plot(square);title('50个随机序列的均方值');四、实验结果及分析由上结果可知：将图中的计算结果平均后，得到的结果为：产生的50个点的随机序列均值的平均值为：EX=0.0090197；产生的50个点的随机序列方差的平均值为DX=1.0078；产生的50个点的随机序列均方值的平均值为RMS=1.0087。
Levinson-Durbin-算法-和-伯格(Burg)算法
~数字信号仿真实现—题目: Levinson-Durbin 算法和伯格(Burg)算法讲课老师:学生姓名：\所属院系：信息科学与工程学院专业：信息与通信工程学号：·完成日期：2015/5/12Levinson-Durbin 算法%检验Levinson-Durbin算法(clear;%清除内存变量clc;%清屏close;%===========================================================%估计2阶自回归模型的功率谱！%步骤1：建立信号模型，产生观测数据xvar = 1;u = var*randn(1,1000);%产生均值为0，方差为1的高斯白噪声u，数据长度为1000 % figure(3);% plot(u);·a0 = [1 0 ];x = filter(1,a0,u);%信号模型，白噪声通过线性系统H(z)=1/(1+a1*z^(-1)+a2*z^(-2))产生信号向量%画出信号x的功率谱,Sxx(exp(j*w))=var/abs(1+sum(ak*exp(-j*w*k)))^2,求和范围为从1到2，系数ak为给定%的模型参量0，w = linspace(-pi,pi,2000);%将-pi到pi均分为2000等分\for mm = 1:2000c = w(mm);S(mm) = var/(abs(1+a0(2:3)*exp(-j*c*(1:2))'))^2;end/subplot(211);plot(w,S,'b');%根据已知参量画出信号x的功率谱xlabel('角频率/rad');ylabel('x的功率谱');title('2阶自回归AR模型的功率谱');\%===========================================================%估计2阶滑动平均模型的功率谱%步骤1：建立信号模型，产生观测数据x·% clear;var1 = 1;u1 = var1*randn(1,1000);%产生均值为0，方差为1的高斯白噪声u，数据长度为1000b0 = [1 1 1];x1 = filter(b0,1,u1);%信号模型，白噪声通过MA(2)阶系统x(n)=u(n)+u(n-1)+u(n-2) :%根据已知信号参量画出信号x的功率谱，S = var*abs(sum(bk*exp(-j*w*k)))^2 for mm = 1:2000c = w(mm);S1(mm) = var1*(abs(b0*exp(-j*c*(0:2))'))^2;end￥subplot(212);plot(w,S1,'b');%根据已知参量画出信号x的功率谱xlabel('角频率/rad');ylabel('x的功率谱');…title('2阶滑动平均MA模型的功率谱');运行结果如下图1 所示：图 1 Levinson-Durbin算法结果图^伯格(Burg)算法clear#%取样点%定义常数值N=32;a(1)=;d2=;f1=;f2=;f3=;*ur=*d2.*randn(1,N);ui=*d2.*randn(1,N);u=ur+ui*i;%定义32个复数点z(1)=u(1);，x(1)=6+z(1);for n=2:Nz(n)=-a(1)*z(n-1)+u(n);x(n)=2*cos(2*pi*f1*(n-1))+2*cos(2*pi*f2*(n-1))+2*cos(2*pi*f3*(n-1))+z(n); end【%定义f范围fmin=;fstep=;fmax=;f=fmin:fstep:fmax;nf=(fmax-fmin)/fstep;t=sqrt(-1);\%初值rxx=0;p0=zeros(1,11);ef=zeros(11,N);eb=zeros(11,N);a=zeros(10,10);for n=1:Nrxx=rxx+(abs(x(n)))^2;end（rxx=(1/N)*rxx;p0(1)=rxx;ef(1,:)=x;eb(1,:)=x;ef(1,1)=0;eb(1,32)=0;%算法p=10;kk=zeros(1,10);for k=1:p]e1=0;e2=0;for n=(k+1):Ne1=e1+ef(k,n)*(conj(eb(k,n-1)));e2=e2+(abs(ef(k,n))^2+abs(eb(k,n-1))^2);kk(k)=(-2)*e1/e2;~ef(k+1,n)=ef(k,n)+kk(k)*eb(k,n-1);eb(k+1,n)=eb(k,n-1)+conj(kk(k))*ef(k,n);endfor i=1:(k-1)a(k,i)=a(k-1,i)+kk(k)*conj(a(k-1,k-i));enda(k,k)=kk(k);p0(k+1)=(1-abs(kk(k))^2)*p0(k);end%功率谱for j=1:nf+1pxx=0;for k=1:ppxx=pxx+a(10,k)*exp(-t*2*pi*f(j)*k);endpxxf(j)=d2/[(abs(1+pxx))^2];pxxf(j)=10*log10(pxxf(j));endplot(f,pxxf)title('伯格(Burg)算法');运行结果如下图2 所示：图 2 Burg 算法结果图。
二维burgers方程的数值求解
二维burgers方程的数值求解二维 Burgers 方程是求解二维水流物理场的实际应用问题，是非线性偏微分方程。
它可以描述例如洋流、湍流或海浪等非稳定状态的特征。
其数值求解方法有：1、采用传统有限差分法：（1）用有限差分法离散二维Burgers方程：采用差分格式将解析形式的二维 Burgers 方程精确地转化成离散的形式，也就是将区域被划分成若干小的网格，然后求解每个网格区域的值。
（2）根据所采取的差分格式，将有限差分方程按照某一特定数值算法求解：不同的算法实现有不同的发展步骤，可以根据算法的具体特点选择数值求解步骤，利用具体的二维 Burgers 方程来求解每个微分方程的方法，目的是为了找出最精确的求解方案，这样可以节省更多的计算资源，并在有限的计算时间内获得更准确的数值求解结果。
2、采用隐式（implicit）有限差分法：（1）采用隐式有限差分法对二维 Burgers方程进行离散：它和传统有限差分法一样，也是先将求解区域分割成若干小区域，然后再利用有限差分格式将解析形式的二维 Burgers 方程精确地转化成离散的形式。
不同的是，这种方法不会受到传统的有限差分法求解的约束，可以解决一定规模问题的更大精度问题，同时消耗更少的计算资源。
（2）采用特定数值算法求解有限差分方程：采用隐式有限差分法和特定数值算法，可以利用全局最优解法以较少的计算时间精确地求解二维 Burgers 方程，消耗更少的计算资源，并且可以在较短的时间内获得准确的求解结果。
3、采用积分方法：（1）采用隐式积分法求解二维 Burgers 方程：采用隐式积分方法求解二维 Burgers 方程时，需要先离散数值化低阶微分形式的微分方程，然后用积分的方法来求解最终的微分方程。
（2）利用特定的数值算法求解积分方程：由于改变网格的精度和求解方法的关系，积分方法求解二维 Burgers 方程有不同的具体发展步骤，利用具体的二维 Burgers 方程来求解每个微分方程的方法可以获得最准确的求解结果。
噪声中正弦信号的现代法频谱分析
实验报告一、实验名称噪声中正弦信号的现代法频谱分析二、实验目的通过对噪声中正弦信号的现代法频谱分析，来理解和掌握现代谱估计的基本概念，以及学会应用现代谱估计以及改进后的方法。
三、基本原理1．参数模型法是现代谱估计的主要内容，思路如下：① 假定所研究的过程)(n x 是由一个输入序列)(n u 激励一个线性系统)(z H 的输出； ② 由已知的)(n x ，或其自相关函数)(m r x ； ③ 由)(z H 的参数来估计)(n x 的功率谱。
2．自回归模型，简称AR 模型，它是一个全极点的模型。
“自回归”的含义是：该模型现在的输出是现在的输入和过去p 个输出的加权和。
此模型可以表现为以下三式：① ∑=+--=pk k n u k n x a n x 1)()()(；② ∑=-+==pk kk z a z A z H 111)(1)(；③ 2121)(∑=-+=pk jwkk jwx e a e P σ。
3．AR 模型的正则方程建立了参数k a 和)(n x 的自相关函数的关系，公式如下：=)(m r x ∑=--pk x k k m r a 1)( 1≥m 时，=)(m r x 21)(σ+-∑=k r a pk x k 0=m 时。
4．Levinson-Durbin 算法：从低阶开始递推，直到阶次p ，给出了在每一个阶次时的所有参数。
公式如下：① 1111/])()()([--=-∑+--=m m k x x m m m r k m r k a k ρ；② )()()(11k m a k k a k a m m m m -+=--；③ ]1[21mm m k -=-ρρ。
5. 自相关法：使用线性预测的方法来计算不同阶数下的预测器系数，使用前向线性预测，预测误差为∑=-+=-=pk k f k n x a n x n x n x n e1)()()(ˆ)()(ˆ，预测均方误差为∑-+==1022}{1)]([p N n ffn eNn e E ，使前向预测误差功率相对AR 参数k a 最小，将反射系数代入Levinson-Durbin 算法即可求解。
张 宏-地震谱分解算法对比与局限性分析
第30卷第6期2007年12月勘探地球物理进展Progress in Exploration G eophysicsVol.30,No.6Dec.,2007收稿日期:2007206225;改回日期:2007209229。
作者简介:张宏(1967—),男,高级工程师,中国地质大学(北京)在读博士,主要从事地震地质资料综合应用和研究工作。
文章编号:167128585(2007)0620409206地震谱分解算法对比与局限性分析张 宏1,2(1.中国地质大学(北京)能源学院,北京100083;2.中国石油化工股份有限公司河南油田分公司勘探开发研究院,河南南阳473132)摘要:地震反射是典型的非平稳信号,地震谱分解实质上就是连续时频谱分析。
谱分解为一非惟一的过程,可供选用的算法很多。
每种算法既有优点,又有缺陷。
为了达到精确应用的目的,选用适当的算法十分重要。
在对目前常用的几种谱分解算法介绍的基础上,结合其中一些方法的谱分解结果,对各种算法的时间、频率分辨率进行了对比分析,指出了各类算法的局限性,并针对不同用途指出了适用的谱分解方法。
关键词:谱分解;时频分析;地震属性;分辨率;储层预测;油气检测中图分类号:P631.4文献标识码:A 传统的信号分析建立在傅里叶变换的基础上。
由于傅里叶分析使用的是一种全局的变换,要么完全在时域,要么完全在频域,因此无法表征信号的时频局部特性,更无法表征信号的瞬时频率特征。
而对于象地震反射信号这样的非平稳信号,更需关注的恰恰就是信号的局部(或瞬时)时频特征。
为此,在傅里叶分析的理论基础上,提出了一系列分析和处理非平稳信号的时频分析技术。
时频分析很早就是地震数据分析和资料处理的重要工具(如频谱分析、频率滤波和子波特征描述等)。
1982年Morlet等[1]首次把短时傅里叶频谱用于地震解释,1999年Partyka等[2]借助先进的计算机技术把时频分析技术转化为一种实用、便捷的解释工具,即所谓的地震“谱分解”[3](Special Decompo sition)技术。
现代数字信号处理题解
现代数字信号处理题解现代数字信号处理复习题⼀、填空题1、平稳随机信号是指：概率分布不随时间推移⽽变化的随机信号，也就是说，平稳随机信号的统计特性与起始时间⽆关，只与时间间隔有关。
判断随机信号是否⼴义平稳的三个条件是：（1）x(t)的均值为与时间⽆关的常数：C t m x =)( （C 为常数）；（2）x(t)的⾃相关函数与起始时间⽆关，即：)(),(),(ττx i i x j i x R t t R t t R =+=；（3）信号的瞬时功率有限，即：∞<=)0(x x R D 。
⾼斯⽩噪声信号是指：噪声的概率密度函数满⾜正态分布统计特性，同时其功率谱密度函数是常数的⼀类噪声信号。
信号的遍历性是指：从随机过程中得到的任⼀样本函数,好象经历了随机过程的所有可能状态,因此,⽤⼀个样本函数的时间平均就可以代替它的集合平均。
⼴义遍历信号x(n)的时间均值的定义为：，其时间⾃相关函数的定义为：。
2、连续随机信号f(t)在区间上的能量E 定义为：其功率P 定义为：离散随机信号f(n)在区间上的能量E 定义为：其功率P 定义为：注意:(1)如果信号的能量03、因果系统是指：对于线性时不变系统，如果它在任意时刻的输出只取决于现在时刻和过去时刻的输⼊，⽽与将来时刻的输⼊⽆关，则该系统称为因果系统。
4、对平稳随机信号，其⾃相关函数为)(τx R ，⾃协⽅差函数为)(τx C ，（1）当0→τ时，有：)(τx R =x D ，)(τx C ＝2x σ。
（2）当∞→τ时，有：)(τx R =2x m ，)(τx C ＝0。
5、⾼斯－马尔可夫随机信号的⾃相关函数的⼀般表达式可表⽰为:||)(τβητ-e R x ＝。
6、⾼斯–马尔可夫信号)(t x 的⾃相关函数为||410)(ττ-e R x ＝，其均值 0)(=∞=x x R m ，均⽅值10)0(==x x R D ,⽅差102==x D σ。
其⼀阶概率密度函数的表达式为：?-=102exp 201)(2x x p π。
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
k 1
( n 1) k e
k 1
k 1
( n 1)
e k ( n) k e
解出
k
( n 1) e
k 1
( n 1)代入之
n k N 1
2 ek 1 ( n) ek 1 ( n 1) [ek 1 ( n) 2 ek 1 ( n 1) 2 ]
R
(2) 迭代从k=1开始
N 1
k 1 n N 1
由：
k
n p N 1 n p
2 ek 1 ( n) ek 1 ( n 1) [ek 1 ( n) 2 ek 1 ( n 1) 2 ]
解出k
(3)由： Ak ( z ) Ak 1 ( z ) k z A
谱线分裂 噪声对AR谱估计的影响
补救方法
(1)采用ARMA谱估计方法 (2)对数据进行滤波,减小噪声
(3)采用高阶AR估计
(4)补偿自相关函数或反射系数估计中噪声 的影响
Burg 算法
引言
Levinson递推的缺点
计反射系数，利用Levinson 递推由反射系数来估计AR参数。使用前后向预 测误差功率估计的平均值最小为最佳准则确定 反射系数，避免了自相关函数的估计。 • 功率估计时用时间平均代替集平均。
定义前向、后向预测值误差功率
N 1
n k
解出p后，再利用Levison递推，由k-1阶 AR 参数计算出k阶AR 参数，具体计算步骤：
e0 ( n ) e0 ( n ) x ( n ) (1)确定初始值：
1 2 0 N
n0
N 1
x 2 ( n)
A0 ( z ) A 0 ( z ) 1
求出e ( n), e ( n)
(5)求确定k阶均方误2k
(6)回到(2)作下一次迭代
2 k
k n N 1
(1 k )
2
2 k 1
不必计算自相关函数
AR估计的异常现象及补救措施
AR估计的异常现象
虚假谱峰 谱线分裂 噪声对AR谱估计的影响 虚假谱峰产生的原因由于估计误差引起ap,i i>p时ap,i不为零 补救方法：降低模型的阶次。不超过N/2
1
R
k 1
( z ) 求出Ak (z
i 1,2,...k
或
a k , i a k 1, i k a k , k i
e k ( n ) 1 (4)由： e k ( n) k
k k
k e k 1 ( n) 1 e k 1 ( n 1)
[e k ( n ) 2 e k ( n ) 2 ] N 1 n k
欲使 min
N 1 2 k n k
将e k ( n) e
求 k 使 min 令 0 k ek ( n) ek ( n) [e k ( n ) e p ( n) ] 0 k k