功-率-谱-估-计要点分析.ppt

为何讨论AR模型？ M=15; %阶数 M1=20; %阶数 %求自相关函数 for m=0:M-1, % 对r赋初值 temp=0; for n=1:N-M, temp=x(n)*x(n+m)+temp; end r(m+1)=(1/N)*temp; end %%%%%%%%%%%%%根据数学模型计算功率谱 q=1; for f=0:0.001:0.5-0.001; w=2*pi*f; t=0; for i=1:p-1 t=t+a(p-1,i)*exp((-1)*j*w*i); end pxx(q)=(1/abs(1+t))^2; q=q+1; end plot(f/(2*pi),pxx,r) f=0.001:0.001:0.5; subplot(2,1,1); axis([0 0.8 -20 120]); plot(f,10*log(pxx)); title(‘自编程序’)；axis([0 0.5 -20 100]) nfft=100; subplot(2,1,2); [pxx1,f1]=pyulear(x,M,nfft,fs);% MATLB自带的求解函数 axis([0 0.4999 -20 120]); plot(f1,10*log(pxx1)); title(用MATLB函数) axis([0 0.5 -20 100])； m=k时，有： 算法的程序编写流程图参见讲义P79或丁P154 图 3.3.1 利用列文森递推法计算功率谱的流程图 表明： 系统函数所有极点均在单位圆内 讨论： fs=1; N=10000; t=0:1/fs:N; A=20^0.5; B=A; C=1;　 x1=A*sin(2*0.15*pi*t); x2=B*sin(2*0.2*pi*t); noise=C*randn(size(t)); x=x1+x2+noise; %被分析信号 练习：自编程序，请理解、消化 a(m+1,m+1)=0; a(1,1)=-r(2)/r(1); xgm(1)=(1-(a(1,1))^2)*r(1); % p=2; while pM p=p+1; temp=0; for i=1:p-2 temp=a(p-2,i)*r(p-i)+temp; end a(p-1,p-1)=-(r(p)+temp)/xgm(p-2); for k=1:p-2 a(p-1,k)=a(p-2,k)+a(p-1,p-1)*a(p-2,p-k-1); end xgm(p-1)=(1-(a(p-1,p-1))^2)*xgm(p-2); end 式中 -(M-1)≤n≤M-1 估计方法： 那么 周期图的窗函数法就是前面BT法的加权协方差谱估计。 又 偏移分析： 估计效果分析： 可得 周期图的窗函数法仍然是有偏估计, 其偏移和wB(m)、w(m)两个窗函数有关。 如果w(m)窗的宽度比较窄，M比N小得多，这样|m|N，则wB(m)~1, 由于w(m) 比wB(m)窄，W(ejw) 的主瓣比WB(ejw)宽，故可以利用窗函数法进一步平滑周期图，减小估计方差；但相应的会增加偏移，降低频率分辨率。 图 4.2.4 利用数据窗减少窄带过程周期图的旁瓣（丁） (a) 没有数据窗； (b) 加哈明数据窗 周期图和谱窗函数卷积得到功率谱，等效于在频域对周期图进行修正，使周期图通过一个线性非频变系统，滤除掉周期图中的快变成分，谱窗函数需具有低通特性 3. 修正的周期图求平均法（Welch法） 主要思想：对Bartlett法进行修正，使之更适合FFT计算。 选择适当的窗函数w(n)，并在周期图计算前直接加进去； 在分段时，可使各段之间有重迭，这样将会使方差减小。 估计方法： 首先把数据长度为N的信号x(n)分成L段，每一段数据长度为M，N=LM； 然后把窗函数w(n)加到每一个数据段上， 求出每一段的周期图，形成修正的周期图； 再对每一个修正的周期图进行平均。 第i段的修正周期图为 i=1, 2, 3, …，L 式中 同样, 将每一段的修正的周期图之间近似看成互不相关, 最后功率谱估计为 对上式求统计平均, 得到 式中 估计效果分析： 估计是渐近无偏的； 这种方法对窗函数没有限制，无论什么样的窗函数均可使谱估计非负； 分段时，为了减少因分段数增加给分辨率带来的影响，可使各段间有重叠，例如重叠50%。 结论：传统的功率谱估计方法无论采取哪一种改进方法，总是以减少分辨率为代价，换取估计方差的减少，提高分辨率的问题无法根本解决。 3.3