自己编写算法的功谱密度的三种matlab实现方法4.docVIP

功率谱密度的三种matlab实现方法 一：实验目的：(1)掌握算法的概念、应用及特点； (2)了解谱估计在信号分析中的作用； (3) 能够利用法对信号作谱估计，对信号的特点加以分析。 的估计的抽样. 认为随机序列是广义平稳且各态遍历的，可以用其一个样本x(n)中的一段来估计该随机序列的功率谱。这当然必然带来误差。由于对采用DFT，就默认在时域是周期的，以及在频域是周期的。这种方法把随机序列样本x(n)看成是截得一段的周期延拓，这也就是周期图法这个名字的来历。 相关法（间接法): 这种方法以相关函数为媒介来计算功率谱，所以又叫间接法。这种方法的具体步骤是： 第一步：从无限长随机序列x(n)中截取长度N的有限长序列列 第二步：由N长序列求（2M-1）点的自相关函数序列。 (2-1) 这里，m=-(M-1)…,-1,0,1…,M-1，MN，是双边序列，但是由自相关函数的偶对称性式，只要求出m=0,。。。,M-1的傅里叶变换，另一半也就知道了。 第三步：由相关函数的傅式变换求功率谱。即 以上过程中经历了两次截断，一次是将x(n)截成N长，称为加数据窗，一次是将x(n)截成（2M-1）长，称为加延迟窗。因此所得的功率谱仅是近似值，也叫谱估计，式中的代表估值。一般取MN，因为只有当M较小时，序列傅式变换的点数才较小，功率谱的计算量才不至于大到难以实现，而且谱估计质量也较好。因此，在FFT问世之前，相关法是最常用的谱估计方法。 三：Burg法： AR模型功率谱估计又称为自回归模型，它是一个全极点的模型，要利用AR模型进行功率谱估计须通过levinson_dubin递推算法由 Yule-Walker方程求得AR的参数:σα1α2…αp。 nfft=10000; %2^n n=0:Fs; x=sin(2*pi*0.2*n)+sqrt(2)*sin(2*pi*0.213*n)+randn(size(n)); X=fft(x,nfft); Pxx=abs(X).^2/length(n); %求解PSD t=0:round(nfft/2-1); f=t/nfft; P=10*log10(Pxx(t+1)); %纵坐标的单位为dB plot(f,P); grid on nfft=200 nfft=1024 nfft=10000 相关法： clear; Fs=1000; %采样频率 n=0:Fs;%产生含有噪声的序列 nfft=1024; xn=sin(2*pi*0.2*n)+sqrt(2)*sin(2*pi*0.213*n)+randn(size(n)); cxn=xcorr(xn,unbiased); %计算序列的自相关函数 CXk=fft(cxn,nfft); %求出功率谱密度 Pxx=abs(CXk); index=0:round(nfft/2-1); f=index/nfft; plot_Pxx=10*log10(Pxx(index+1)); plot(f,plot_Pxx); xlabel(频率); ylabel(功率/DB); grid on; nfft=256 nfft=1024 Burg法： clear Fs=1000 %设置关键变量，可通过调节这些变量观察不同效果 f1=0.2; f2=0.213; nfft=128; %取样点数 p=50; %阶数p应该选择在N/3pN/2 delta=1; m=sqrt(-1); f=0:1/1000:0.5; n=1:Fs; xn=sin(2*pi*f1*n)+sqrt(2)*sin(2*pi*f2*n)+randn(size(n)); figure; plot(n,xn); title(burg时域); xn= xn(:); N=length(xn); ef = xn; eb = xn; a = 1; for l=1:p efp = ef(2:end);%m-1阶前向预测误差 ebp = eb(1:end-1);%m-1阶后向预测误差 num = -2.*ebp*efp;%1km分子多项式 den = efp*efp+ebp*ebp;%1km的分母多项式 k(l) = num ./ den;%计算反射系数 % 更新前向和后向预测误差 ef = efp + k(l)*ebp;%各阶前向预测误差 eb = ebp + k(l)*efp;%各阶后向预测误差 % 计算模型参数 a=[a;0] + k(l)*[0;c