kuirs_doa_multiph.m

ofdm多径仿真各径的到达角度DOA 有详尽的注释 欢迎交流 [e-mail]

%-------------------------------------------------------------------
%                OFDM TOA DOA ESTIMATION
%
%                   KUIRS CHEN 2009.02%
%                           CUMT
%                     [e-mail]
%--------------------------------------------------------------------
clear all; clc;
fftsize=128; cp=fftsize/8;         % FFT 长度 循环前缀长度
lamda=3e8/2.4e9;d=lamda*0.5;         % 载波波长 天线单元间距
%============================ ======================================
M =4; data = randint(fftsize,1,M); % 调制阶数,随机整数Message signal
mod_data = qammod(data,M,0);       % QAM调试
% scatterplot(mod_data);           % 散点图    ;
train=creat_cazac(fftsize,127);    %  训练数列 CAZAC序列
txd=[train mod_data];              % 插入训练序列
ofdm=ifft(txd,fftsize);            % OFDM 调制
txx=[ofdm((fftsize-cp+1):fftsize,:);ofdm]; % 插入循环前缀
tx=[txx(:,1); txx(:,2)];                   % 并－－－－》串
%=================================================================
an=8;                                     % 天线数量，决定峰值的尖锐程度
ch2=[0.8441 zeros(1,2) 0.5 0 0 0.20 0 0 0.16661]; % 信道脉冲响应
%=================================================================
theta=d2r([60 20 22 59]);                          % 各径到达的角度
a_t=exp(j*2*pi/lamda*(1:an)'*d*cos(theta));      % 导向矩阵 an×theta
a_m=[a_t(:,1) zeros(an,2) a_t(:,2)...
    zeros(an,2) a_t(:,3) zeros(an,2) a_t(:,4)];  % 导向矩阵多径处理

snr=30;
%  SNR=[4:2:16];
%  for s =1:length(SNR)
%      an=SNR(s);
%===============================================================================   
ch_m=[];
for i=1:an
    ch_m=[ch_m;ch2]; % 天线接收脉冲响应矩阵
end
chh=a_m.*ch_m;       % 信道脉冲响应矩阵乘上导向矩阵
%================================================================================
% hegecishu=0;doaerr=0;simnum=100; % 实验次数
% for s=1:simnum




for sample=1:100  % 快拍数
    for i=1:an
        rdata(i,:) = filter(chh(i,:),1,tx);          % 信号通过各天线
        rdata(i,:) = awgn(rdata(i,:),snr,'measured');  % 加入AWGN噪声
        TE = rdata(i,17:144);     % 去除循环前缀，取出训练序列的接收量 
        p(i,:)=fft(TE,128);       % 结合信号做FFT变换，到频域  
    end
    sa_m(:,:,sample)=p;
end
sam=mean(sa_m,3);       % 对OFDM个载波数据进行平均
%==================== MUSIC  ALGORITHM ==============================================


RCC=sam*sam'/length(sam);                   % 接收信道矩阵相关 修改后为（AN×AN）矩阵
flag=2;                         %  1：正常的方法 2：改进的方法
JJ=fliplr(eye(length(RCC)));    % 修正矩阵
MR=(RCC+JJ*RCC.'*JJ)/2;         % 修正后的相关矩阵

if flag==1
    [U,V,uu]=svd(RCC);
else
    [U,V,uu]=svd(MR);
end

Un=U(:,[5:end]);
dtheta=0.005; % doa 搜索步长
ww=[0 100];     % MUSIC搜索范围，减少仿真计算量
n=[ww(1):dtheta:ww(2)-dtheta];
for i=1:length(n)
    aa=zeros(an,1);
    for m=1:an
        aa(m,1)=exp(j*(m-1)*pi*cos(n(i)*pi/180));
    end
    P(i)=-10*log(abs(aa'*Un*Un'*aa));
end
P=P./sqrt(sum(P.^2));


%============ 搜索DOA峰值  ===============================================
fanwei=[1: length(P)];
[peak_h,index_h]=findpeak(P(fanwei)); % 搜索范围 缩小范围减少计算量
t_h=(index_h-1)*dtheta;
Lp=4;   %多径的数量，即有多少个不同的入射角度
[DOAa,Atten]=selectpath(peak_h,t_h,Lp);
% ==========================================================================
disp(['天线数量：' num2str(an) ' 信噪比：' num2str(snr)]);
disp(['估计的DOA：' num2str(DOAa)]);

% s
% temp=min(DOAa)-10;
% if abs(temp)<10
%     doaerr=(temp-20)^2+doaerr;
%     hegecishu=hegecishu+1
% end
% end
% 
% rmse=sqrt(doaerr/hegecishu);
% disp(['DOA的RMSE误差：' num2str(rmse)]);


% 
i=1:length(n);
plot(n(i),P(i));hold on;
grid on; title(['MUSIC DOA (SNR= ' num2str(snr) ' dB )']);
xlabel('来波方向角');ylabel('P(dB)')



% DOA_M(s,:)=DOA;
%  end
% DOA_M=sort(DOA_M,2)
% % plot(SNR,sqrt((x.^2)))