利用Pisarenko谐波分解恢复理论推导谐波过程matlab实现

function fangzhen1_ls() pe=6; %生成增广矩阵Re M=100; loops=20; %pe=25; %Re=zeros(M,pe+1); fvMatrix=zeros(pe/2,loops); arMatrix=zeros(pe,loops); for loop=1:1:loops % w = wgn(1,2000,0); w = randn(2000,1); for n=1:1:128 x(n)=(20^(1/2))*sin(2*pi*0.2*n)+(2^(1/2))*sin(2*pi*0.213*n)+w(n+50*loop); %x(n)=(20^(1/2))*sin(2*pi*0.2*n); end Rxx=xcorr(x,'unbiased'); for i=1:1:M for j=1:1:pe+1 Re(i,j)=Rxx(pe+i+1-j+128); end end kk=pe; %最小二乘法估计ARMA模型的参数 b=-1*Re(:,[1]); A=Re(:,[2:1:kk+1]); %ar_lsline=(A'*A)\(A'*b); ar_lsline=inv(A'*A)*A'*b; arMatrix(:,loop)=ar_lsline; %求谐波频率 for j=2:1:kk+1 fc(j) = ar_lsline(j-1); end fc(1) = 1; fz=roots(fc); count=1; for j=1:2:kk fw(count)=atan(imag(fz(j))/real((fz(j)))); fv(count)=fw(count)/(2*pi); count=count+1; end fv=abs(fv); fv=sort(fv); fvMatrix(:,loop)=fv'; end %分别求AR参数估计和频率估计的均值和方差 arMatrix ar_mean = mean(arMatrix')' ar_var = var(arMatrix')' fvMatrix fv_mean = mean(fvMatrix')' fv_var = var(fvMatrix')'