频率估计matlab

%频率估计—music %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %------------------------------------------------------------------------ %初始条件 %------------------------------------------------------------------------ clear; clc; for L=1:20; n=1:128; xn=2*cos(2*pi*0.05*n)+3*cos(2*pi*0.40*n)+1.2*cos(2*pi*0.42*n)+sqrt(0.32)*randn(size(n)); p=6; %阶数； m=p*2; %------------------------------------------------------------------------- %1、由xn求r(0)—r(m-1)，构成自相关矩阵Rxx %------------------------------------------------------------------------- N=length(xn); rxx1=xcorr(xn,'biased'); %求自相关函数； rxx=rxx1(N:m+N-1); %取r(0)—r(m-1)； Rxx=toeplitz(rxx); %m阶自相关矩阵 %------------------------------------------------------------------------- %2、对Rxx特征分解，求出主、次特征值，并构成U1和U2 %------------------------------------------------------------------------- [V D]=eig(Rxx); ev=eig(Rxx); %矩阵分解； [s i]=min(ev); U2=V(:,1:m-p); U1=V(:,m-p+1:m); %------------------------------------------------------------------------- % 3、计算Smusic(f)=1/(e’*U2*U2’*e) %------------------------------------------------------------------------- f=0.0; %f=0:0.0001:0.5; so k is range from 1 to 5000; df=0.0001; df=0.001; %df 画图精度 for k=1:0.5/df; %此循环 将Smusic的值算出来 for i=1:m; %定义e(i) e(i)=exp(-2*pi*(i-1)*f*j); end Smusic(k)=1/abs(e*U2*U2'*e'); f=f+df; end %------------------------------------------------------------------------- %4、求S取最大值的p个峰值对应的横坐标fx，即k*df,求k值; %------------------------------------------------------------------------- f1=0.1/df; f2=0.41/df; f3=0.5/df; s1=max(Smusic(1:f1)); s2=max(Smusic(f1:f2)); %There is a warning here,but I don't know how to correct it. s3=max(Smusic(f2:f3)); %A warning here,too. k1=find(Smusic==s1); k2=find(Smusic==s2); k3=find(Smusic==s3); fx(L,1:3)=df*[k1 k2 k3]; %------------------------------------------------------------------------- %5、画图20合一； %------------------------------------------------------------------------- f=0.0:df:0.5-df; plot(f,Smusic); xlabel('f/Hz'); ylabel('频谱'); title('Music方法-20次循环叠加'); hold on end clc; %------------------------------------------------------------------------- %end %------------------------------------------------------------------------- %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%