oja_mpca_music.m

用基于Oja准则的PCA神经网络方法实现MUSIC算法

close all
clc
clear

c=3e8;   %光速
f0=1e8;  %中心频率100兆
%BW=2e7;                                                                %带宽为20兆
bi=0.45;                                                                %阵间距和波长的比
num_sen=16;                                                             %阵元数
snr=[-10,-10,-10];                                                         %信噪比 db
seta=[-10,0,20].*pi/180';
num_it=40;
fs=4e8;                                                                  %采样率为400M
Ts=1/fs;
degrad=pi/180;                                                           %将角度制化为弧度制 
num_snap=256;
num_p=3;
sigma=0.2;
N=sigma*(randn(num_sen,num_snap)+j*randn(num_sen,num_snap));
snr=sqrt(2*sigma.^2*10.^(snr/10));
t=(1:num_snap)*Ts;

s1=snr(1).*exp(j*2*pi*f0*t+pi/4);
s2=snr(2).*exp(j*2*pi*0.97*f0*t+pi/4);
s3=snr(3).*exp(j*2*pi*0.94*f0*t+pi/4);
S=[s1;s2;s3];
A1=zeros(num_sen,num_p);
A1=exp((-j*2*bi*pi*(0:num_sen-1)).'*sin(seta));
X=A1*S+N;
%  eX=mean(X.');
Rxx=X*(X/num_snap)';
%  for ii=1:num_snap
%      X(:,ii)=X(:,ii)-eX.';
%  end
% eX1=mean(X.');
[E1,T1,V1]=svd(Rxx); 
[E2,T2]=eig(Rxx);

lm=norm(X);                                          
eta=1/(1.2*lm.^2);   

%eta=0.75;
%%%%%%%%学习速率?????
%W1=rand(num_sen,num_sen)+j*rand(num_sen,num_sen);
W1=0.1*rand(num_sen,num_sen);                                               %%%%%%%%%%%%%%%%%%%%%%%付初值%%%%%%%%
W2=zeros(num_sen,num_sen);
Y=zeros(num_sen,num_sen);
CC=zeros(num_sen,num_sen);

E=zeros(num_it,1);

for n=1:num_it
    
    %eta=(1-n/(7.2*num_it))*eta;
    Y=W1*Rxx;
 for ii=1:num_sen
     CC=zeros(num_sen,num_sen);
    for jj=1:num_sen
        CC=CC+W1(jj,:)'*Y(jj,:); 
    end
    BB = Rxx - CC;
    W2(ii,:)=W1(ii,:)+eta.*Y(ii,:)*BB';
 end
  EX(n)=(norm(W2-W1));
  
  W1=W2;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
end
W1=W1.';
tempW=W1*W1'; 
tempE=E1*E1';
tempE2=E2*E2';

Vs=W1(:,1:num_p);          
Vu=W1(:,num_p+1:num_sen);                                 %噪声子空间
i=sqrt(-1);
th2=[-90:0.1:90]';                                        %线性信号角度范围
A3=exp(2*bi*j*pi*(0:num_sen-1).'*sin(th2'.*pi/180));      %[-90 90]范围内的相位差组成的向量
num=diag(A3'*A3);      
Ena=Vu'*A3;
%den=diag(Ena'*Ena); 
den=diag(A3'*(eye(num_sen,num_sen)-Vs*Vs')*A3); 
doa=num./den;                                             %MUSCI估计值
doa1=10*log(doa);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%做出空间谱图
%semilogy(th2,doa); 
figure(1)
ZX=1:num_it;
plot(ZX,EX);
title('Oja权值W收敛图');    %图表标题，x轴y轴标签及范围
xlabel('   迭代次数   ');
%xlabel('   Iteration    ');
ylabel('误差');
%xlabel('   Angle   ');
%ylabel('Spatial Spectrum(dB)');
figure(2)

%ylabel('Weight Error');
%legend('SNR=15;  m=8;  N=1024;')
plot(th2,doa1);
grid on;%y轴对数坐标曲线
title('Oja对称子空间方法');    %图表标题，x轴y轴标签及范围
xlabel('   角度   ');
ylabel('空间谱(dB)');
hold on
[E,T,V]=svd(Rxx);                     %计算特征值
                                          % 此处应添加信号源数目估计程序AIC
tt=diag(T);
for k=1:(num_sen);
    he=sum(tt(k+1:num_sen))/(num_sen-k);
    ji=prod(tt(k+1:num_sen))^(1/(num_sen-k));
   % MDL(k)=N*(m-k)*log(he/ji)+0.5*k*(2*m-k)*log(N);
    AIC(k)=2*num_snap*(num_sen-k)*log(he/ji)+2*k*(2*num_sen-k);
end;
[c num_p]=min(AIC);
 
Vs=E(:,1:num_p);          
Vu=E(:,num_p+1:num_sen);                               %噪声子空间
i=sqrt(-1);
th2=[-90:0.1:90]';                                       %线性信号角度范围
A3=exp(-2*bi*j*pi*(0:num_sen-1).'*sin(th2'.*pi/180));  %[-90 90]范围内的相位差组成的向量
num1=diag(A3'*A3);      
%Ena=Vu'*A3;
den1=diag(A3'*(eye(num_sen,num_sen)-Vs*Vs')*A3);  
%en1=diag(Ena'*Ena);     
doa1=num1./den1;                                           %MUSCI估计值
doa1=10*log(doa1);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%做出空间谱图
%semilogy(th2,doa); 
plot(th2,doa1,'-.');
legend('Oja-K','MUSIC')