clc;clear;
%%%%%%%%%%%%%%用Compressive?Sensing方法利用估計MIMO雷達的DOA
M=40;????????????????????????????????????%發射陣元個數
N=5;?????????????????????????????????????%接收陣元個數
L=256;???????????????????????????????????%快拍數目
K=3;?????????????????????????????????????%稀疏度
SNR=0;??????????????????????????????????%信噪比
%f=8.62e12;???????????????????????????????%載波頻率
%c=3e8;???????????????????????????????????%傳播速度
%lamda=c/f;???????????????????????????????%波長
dt=0.5;????????????????????????????%發射陣元間距取半個波長
dr=0.5;????????????????????????????%接收陣元間距取半個波長
%%%%%有3個目標方位2°3°4°
theat0=[-5?1?7];??????????????????????????%目標與陣列的夾角
%d=[1e9?2e8?3e10];????????????????????????%目標距離坐標原點即第一個發射天線的距離
beata=[5?3?2?4];???????????????????????????%目標的反射系數
Q=length（theat0）;
s0=hadamard（L）;??????????????????????????%產生一組正交信號
for?m=1:M
????s（:m）=s0（:m）;??????????????????????%M個發射信號源
end
n=sqrt（10^（-SNR/10）/2）;
for?q=1:Q
????a（:q）=exp（-j*2*pi*dt*sin（theat0（q）*pi/180）*[0:M-1]‘）;%發射陣列引導矢量
????b（:q）=exp（-j*2*pi*dr*sin（theat0（q）*pi/180）*[0:N-1]‘）;%接收陣列引導矢量
end

noise=n*randn（NL）+j*n*randn（NL）;????????%N個接收通道L次快拍的噪聲
r=zeros（NL）;
for?q=1:Q
????r=r+beata（q）*b（:q）*a（:q）.‘*s.‘;
end
r=r+noise;????????????????????????????????%接收的信號
rr=reshape（r.‘N*L1）;
%%%%%CS定理
theat=-10:0.1:10;
P=length（theat）;???????????%將角度空間稀疏化
for?p=1:P
????aa（:p）=exp（-j*2*pi*dt*sin（theat（p）*pi/180）*[0:M-1]‘）;??????%發射陣列引導矢量
????bb（:p）=exp（-j*2*pi*dr*sin（theat（p）*pi/180）*[0:N-1]‘）;??????%接收陣列引導矢量
????psi（:p）=reshape（（bb（:p）*aa（:p）.‘*s.‘）.‘N*L1）;??????????%稀疏域的基矩陣???????????????????
end

KK=60;
%for?n=1:N
phi1=randn（KKN*L）+j*randn（KKN*L）;?????????????????????????????%測量矩陣（高斯分布白噪聲）
phi2=randn（KKN*L）+j*randn（KKN*L）;
phi=（phi1+phi2）/2;
T=phi*psi;??????????????????????????????????????????????????????%恢復矩陣（測量矩陣*基矩陣）
y=phi*rr;???????????????????????????????????????????????????????%獲得線性測量
%end

hat_y=zeros（1P）;??????????????????????????%待重構的譜域（變換域）向量?????????????????????
Aug_t=[];??????????????????????????????????%增量矩陣（初始值為空矩陣）
r_n=y;?????????????????????????????????????%殘差值


for?times=1:K;?????????????????????????????%迭代次數（有噪聲的情況下該迭代次數為K）
????for?col=1:P;???????????????????????????%恢復矩陣的所有列向量
????????product（:col）=abs（T（:col）‘*r_n）;???%恢復矩陣的列向量和殘差的投影系數（內積值）?
????end
????[valpos]=max（product）;????????????????%最大投影系數對應的位置
????Aug_t=[Aug_tT（:pos）];????????????????%矩陣擴充
????T（:pos）=zeros（KK1）;???????????????????%選中的列置零（實質上應該去掉，為了簡單我把它置零）
????aug_y=（Aug_t‘*Aug_t）^（-1）*Aug_t‘*y;????%最小二乘使殘差最小
????r_n=y-Aug_t*aug_y;?????????????????????%殘差
????pos_array（times）=pos;??????????????????%紀錄最大投影系數的位置
end
hat_y（pos_array）=aug_y;????????????????????%重構的譜域向量
hat_x=psi*hat_y.‘;?????????????????????????%重構得到時域信號
orignal=zeros（P1）;
for?q=1:Q
????orignal（theat0（q）*10+101）=beata（q）;
%?????orignal（theat0（q）+1）=beata（q）;
end
p=（-P+1）/2:（P-1）/2;
plot（（p）/10abs（hat_y）‘b-‘）;
%?plot（p-1abs（hat_y）‘b-‘p-1orignal‘r:‘）;
title（‘壓縮感知估計波達角‘）;
legend（‘估

?屬性????????????大小?????日期????時間???名稱
-----------?---------??----------?-----??----
?????文件????????3604??2020-04-09?09:28??4_11_基于壓縮感知理論的DOA估計\doa_cs_zhangjihong.m
?????目錄???????????0??2020-03-21?16:01??4_11_基于壓縮感知理論的DOA估計\