clear;
clc;
close all;
%-------------<<<<参数设置>>>>-------------%
M=10; %%阵元数目
L=1024*2*2*2;
N=3; %%信源数目
ima=sqrt(-1);
thetaa=[-30 -10 -40 ];
snr=10;
C=3e8;
fs=4e8; %%采样率
Ts=1/fs;
%--------------<<<<信号构造>>>--------------%
t=(0:L-1)*Ts;
f=[1e6;1.5e6;2e6;]; %%五个信号频点2.5e6;3e6
%S=randn(3,L); %%五个信号
S=exp(j*2*pi*f*t);
% S=ones(3,L);
f0=100e6; %%中心频率
lamda=C/f0;
d=0.5*lamda;
%%%导向矢量
a=[0:M-1]';
A=exp(j*2*pi*a*d*sin((thetaa)/180*pi)/lamda);
%%%信号产生
X0=A*S;
X=awgn(real(X0),snr)+j*awgn(imag(X0),snr); %%加噪声
Rx=X*X';
S0=[1,zeros(1,9)]';
wopt=inv(Rx)*S0*inv(S0'*inv(Rx)*S0);
X_out=wopt'*X;
figure(1)
plot(((0:L-1)-L/2)*fs/L,20*log10(fftshift(abs(fft(X_out)))));
hold on;
plot(((0:L-1)-L/2)*fs/L,20*log10(fftshift(abs(fft(X(1,:))))),'r');
[V,D]=eig(Rx); %求矩阵的特征值
E=diag(D);
u=0.0002;
w=zeros(9,1);
w_1=[];
for i=1:L
d0=X(1,i);
xn=X(2:10,i);
e=d0-w'*xn;
w=w+u*xn*e';
% w=w-2*u*(eye(10)-S0*S0'/(S0'*S0))*X*X'*w;
w_1=[w_1,w];
end
w_opt=[1;-w];
%w_opt=w;
X_out2=w_opt'*X;
plot(((0:L-1)-L/2)*fs/L,20*log10(fftshift(abs(fft(X_out2)))),'g');
figure(2)
subplot(411)
plot(real(w_1(1,:)))
subplot(412)
plot(imag(w_1(1,:)))
subplot(413)
plot(real(w_1(2,:)))
subplot(414)
plot(imag(w_1(2,:)))
F_mvdr=zeros(1,181);
theta=-90:90;
for i=1:length(theta)
a_theta=exp(1i*2*pi*d/lamda*a*sin((theta(i))/180*pi)); %%
F_mvdr(i)=wopt'*a_theta;
F_mvdr1(i)=w_opt'*a_theta;
end
F_mvdr=20*log10(abs(F_mvdr)/max(F_mvdr));
F_mvdr1=20*log10(abs(F_mvdr1)/max(F_mvdr1));
figure(6)
subplot(211)
plot(theta,F_mvdr);
subplot(212)
plot(theta,F_mvdr1);
% [V1,D1]=eig(Rx1); %求矩阵的特征值
% %提取D矩阵的对角线构成一个一列矩阵
% E1=diag(D1);
% Un=V(:,1:M-N); %噪声子空间
% Un1=V1(:,1:M-N);
% %--------------<<<<谱峰搜索>>>-------------%
% theta=-90:0.2:90;
% for ff=1:length(theta)
% a_theta=exp(1i*2*pi*d/lamda*a*sin((theta(ff)-20)/180*pi)); %%
% Pmusic(ff)=10*log10(1/abs((a_theta'*Un*Un'*a_theta))); %%Music搜索
% Pmusic1(ff)=10*log10(1/abs((a_theta'*Un1*Un1'*a_theta)));
% Capon(ff)=10*log10(1/abs(a_theta'*inv(Rx)*a_theta)); %%Capon算法搜索
% end
%
% figure(1);
% plot(theta,Pmusic);
% xlabel('DOA(°)'),ylabel('Pmmusic(db)');title('经典MUSIC算法');
% grid on
% hold on;
% plot(theta,Pmusic1,'r');
% xlabel('DOA(°)'),ylabel('Pmmusic1(db)');title('经典MUSIC算法');
%
% figure(2);
% plot(theta,Capon);
% xlabel('DOA(°)'),ylabel('谱峰(db)');title('Capon算法')
% grid on
%
% thetab=[60 55 80 82 110 76];
%
% A1=exp(j*2*pi*a*d*sin((90-thetab)/180*pi)/lamda);
% S1=randn(2*N,L);
% %%%信号产生
% X00=A1*S1;
% X11=awgn(real(X00),snr)+j*awgn(imag(X00),snr); %%加噪声
%
%
%
% r=X11*X11'/L;
% R=pinv(r);
%
%
% c=A1;
% f=[1 0 0 0 0 0]';
% % c=[a1theta a4theta a5theta a6theta];
% % f=[1 0 0 0]';
% Wopt=R*c*inv(c'*R*c)*f;
% F_mvdr=zeros(1,181);
% theta=0:180;
% for f=1:length(theta)
% a_theta=exp(1i*2*pi*d/lamda*a*sin((90-theta(i))/180*pi)); %%
% F_mvdr(i)=Wopt'*a_theta;
% end
% F_mvdr=20*log10(abs(F_mvdr)/max(F_mvdr));
%
%
%
% % f0=1e6;
% % f1=2e6;
% % f2=3e6;
% %
% % fs=30e6;
% % ts=1/fs;
% % t=(1:samp)*ts;
% % SIR_1=0;
% % SIR_2=0;
% % SNR1=0;
% %
% % S0=sin(2*pi*f0*t);
% % S1=10^(-SIR_1/10)*sin(2*pi*f1*t);
% % S2=10^(-SIR_2/10)*sin(2*pi*f2*t);
% %
% %
% %
% % S=[S0;S1;S2];
% % noise=10^(-SNR1/10)*(randn(N,samp)+ima*randn(N,samp));
% % X=c1*S+noise;
% %
% %
% % Y=Wopt'*X;
% % figure(1)
% % plot(t,S0);
% % title('S0')
% %
% % figure(2)
% % plot(t,real(X(1,:)));
% % title('X')
% %
% % figure(3)
% % plot(t,real(Y));
% % title('Y')
% %
% %
% %
% % %
% % % Rd=X*S0'/samp; %未知方向时是否可以估计的出来
% % % %Rd=A(:,1); %此法也可以,必须知道方向
% % % R=X*X'/(samp*1);
% % % Wc_MMSE=pinv(R)*Rd;
% % % Wc_MMSE=Wc_MMSE/norm(Wc_MMSE);
% % %
% % %
% % % r1=X*X'/samp;
% % % R1=pinv(r1);
% % % Wopt1=R1*c*inv(c'*R1*c)*f;%%*f
% % %
% % % Y1=Wopt1'*X;
% % %Y1=Wc_MMSE'*X;
% % % figure(4)
% % % plot(t,real(Y1));
% % % title('Y1')
% %
% % figure(5)
% % subplot(311)
% % plot(fs*((1:1024)-1024/2)/1024,20*log10(fftshift(abs(fft(X(1,:),1024)))));
% % title('X')
% % subplot(312)
% % plot(fs*((1:1024)-1024/2)/1024,20*log10(fftshift(abs(fft(Y,1024)))));
% % title('Y')
% % % subplot(313)
% % % plot(fs*((1:1024)-1024/2)/1024,20*log10(fftshift(abs(fft(Y1,1024)))));
% % % title('Y1')
% %
% % % for i=1:length(theta)
% % % a_theta=exp(ima*2*pi*d_lamad*sin(theta(i)*pi/180)*[0:N-1]');
% % % F_mvdr1(i)=Wopt1'*a_theta;
% % % end
% % % F_mvdr1=20*log10(abs(F_mvdr1)/max(F_mvdr1));
% %
% % figure(6)
% % subplot(211)
% % plot(theta,F_mvdr);
% % % subplot(212)
% % % plot(theta,F_mvdr1);
% %
% %
% %
% %
% %
% %
% %
% % % plot(abs(fft(Wopt,1024)));
% % % Wopt=R*c*inv(c'*R*c);
%
%