谐波频率估计的总体最小二乘方法TLS_ESPRIT.rar

%------------TLS_ESPRIT------------------------ clc; clear; number=100; %运行次数 mistake=0; %错误次数 l=0; MSE=0; %均方误差和初值置0 s=9; for num=1:number N=100; % 采样点数 p=4; % 谐波分量个数 m=10; % 选定时间窗宽度 s=m-1; %SNR = 30; % 信噪比 omega1=-0.19*pi; %信号1的角频率 omega2=-0.37*pi; %信号2的角频率 omega3=0.19*pi; %信号3的角频率 omega4=0.37*pi; %信号4的角频率 s1=20; %信号1的幅值 s2=0.2; %信号2的幅值 s3=20; %信号3的幅值 s4=0.2; %信号3的幅值 Ts=0.01; %采样周期 x = zeros(1,N); for k = 0:(N-1) x1=s1*exp(j*omega1*k*Ts); x2=s2*exp(j*omega2*k*Ts); x3=s3*exp(j*omega3*k*Ts); x4=s4*exp(j*omega4*k*Ts); x(k+1)=x1+x2+x3+x4; %x(k+1) = s*exp(j*omega*k*Ts).'; end %y = awgn(x,SNR); % 加性白噪声 x5 =0.5*randn(1,N); %产生均值为0，方差为1的高斯噪声 y=x+x5; r=xcorr(y,m); %计算自相关函数r(0),r(1),...,r(m) Rxx=zeros(m,m); Rxy=zeros(m,m); for k=1:m %生成矩阵Rxx,Rxy Rxx(k,:)=r(m+2-k : 2*m+1-k); Rxy(k,:)=r(m+3-k : 2*(m+1)-k); end d=eig(Rxx); %计算Rxx的特征值 sigma = mean(d(p+1:m)); %对于m>p的特征值取平均作为噪声方差σ^2的估计 for i=1:m % 产生特殊矩阵Z for j=1:m if (i-j==1) Z(i,j)=1; else Z(i,j)=0; end end end Cxx = Rxx - sigma*eye(m); %计算Cxx和Cxy Cxy = Rxy - sigma*Z; [U,D,V]=svd(Cxx); %对Cxx奇异值分解 U1=U(:,[1:s]); V1=U(:,[1:s]); D1=D([1:s],[1:s]); C=U1'*Cxy*V1; gamma = eig(D1,C); %计算矩阵束(D1,C)的广义特征值 gammaeff=0; % 去除广义特征值里偏离较大的点 for i=1:s if (abs(gamma(i))>=0.99 & abs(gamma(i)<=1.01)) gammaeff=[gammaeff gamma(i)]; end end gammaeffective=gammaeff(2:length(gammaeff)); ome=angle(gammaeffective); %将特征值转换成弧度角，即估算出omega omeg=sort(ome); if length(omeg)~=4 mistake=mistake+1; %统计估算出的omega个数不等于给定频率个数的次数 omeg=[0 0 0 0]; %若个数不为给定频率次数，置0 end omegaout(num,:)=omeg(1:p); %将估算出的omega由大到小重新排列并存入矩阵 %-------------------------------------------------------------------------- theta=0:pi/100:2*pi; x1=sin(theta); y1=cos(theta); %绘制单位圆及给定频率 plot(x1,y1); hold on plot(theta,tan(0.19*pi)*theta); hold on plot(theta,tan(0.37*pi)*theta); hold on omega=sort([omega1 omega2 omega3 omega4]); % 给定频率分量并按从小到大排序 j=sqrt(-1); plot(exp(j*omega),'ro') % 给定频率 axis([-1 1 -1 1]); %-------------------------------------------------------------------------- hold on plot(gammaeffective,'x') % 将估计出的频率绘制出来 %hold off end for i=1:number if omegaout(i,:)~=0 l=l+1; omegaou(l,:)=omegaout(i,:); MSE=MSE+(omegaout(i,:)-omega)*(omegaout(i,:)-omega)'; end RMSE=sqrt(MSE/(number-mistake)); %计算number次实验均方误差 end omega omegaou hold on title(strcat('TLS-ESPRIT算法',' 实验次数=',num2str(number),' RMSE=',num2str(RMSE))); xlabel ('红色圈代表给定频率，蓝色×代表估算出的频率'); figure; plot(real(y)); title(strcat('信号波形图',' a1=',num2str(s1/2),' a2=',num2str(s2/2),' 采样周期=',num2str(Ts),' 噪声方差=0.5')); xlabel ('采样点数'); %mistake RMSE