一维线阵干涉仪测角 parfor多线程循环 matlab仿真

duibi表示对比实验，其他部分为信号处理一维线阵的干涉仪测向仿真

clc;
clear all;
close all;

%% 参数
order = 1; %BPSK=1;QPSK=2;
Rs = 10e6; %符号速率
fs = 40e6; % 采样率
osr = fs/Rs;
c = 3e8;
fc = 1.5e9;  %载波频率
N = 3; %阵元个数 除去0阵元
mx = 2;
my = 2;
lambda = c / fc;  %波长
d_each = lambda / 2;  %天线阵元间距
dx =0:d_each:(N-1)*d_each;  %x方向阵元位置
dy =0:d_each:(N-1)*d_each;  %y方向阵元位置
span = 6;
theta = [67];
phi = [56];

EbN0 = -20:10;
%% 发射机
data = randi([0 1],1,512);
data_bio = 1-2*data;
data_up = upsample(data_bio,osr);
hrc =rcosdesign(0.35,span,osr,'sqrt');
baseband = conv(data_up,hrc); %不'same' 滤波器有延时

%% 产生各阵元相位差
for m=1:N
    A(m,1)=exp(-1i*pi*m*( sind(theta) ));  %导向矢量
end
A1 = [1;A]; %方向矩阵
sigdata_in = A1*baseband;

plotData1 = zeros(length(EbN0),200);

%% 信道加噪
for EbN0_index = 1:length(EbN0) 
    
    parfor i_index = 1:200
        %信道 +噪声
        disp([EbN0_index,i_index]);
        [theta_es_out,theta_ess_out]= F_GAN(EbN0,EbN0_index,sigdata_in,order,fs,Rs);
        plotData1(EbN0_index,i_index) = abs(theta_es_out - theta);
        
    end            
end

%% 图
figure;
plot(EbN0,mean(plotData1,2),'b');grid on;

%%
function [theta_es_out,theta_ess_out]= F_GAN(EbN0,EbN0_index,sigdata_in,order,fs,Rs)
    SNR = EbN0(EbN0_index) + 10*log10(2) + 10*log10(order) - 10 * log10(fs/Rs);
    for s = 1:4
        sigdata_in_awg(s,:) = awgn(sigdata_in(s,:),SNR,'measured');  %接收数据
    end

    
    %% 接收
    %% 干涉仪
    
    for ii = 1:3
        k(ii)= sigdata_in_awg(1,:)*sigdata_in_awg(ii+1,:)'; %测角输入数据 %天线(1,2) (1,3) (1,4)  %用降采样之后的或者接收数据没啥影响
    end
    phrase = angle(k); %获取相位差
    
    for m = 1:3
        theta_es(m) = rad2deg(asin( phrase(m)/(pi * m) ));  % 获得估计角度
    end
    
    
    theta_es_out  = theta_es(1);
   
end