%% [预处理]
clc;
clear;
close all;

%% [配置参数]
N=3000000;             %数据点数（个）
SNR_dB=0:1:45;         %信噪比（dB形式）
SNR=10.^(SNR_dB/10);   %信噪比（一般形式，Eb/N0）
Ps=1;                  %信号功率，功率为1时其dB形式为0。
Pn=Ps./SNR;            %噪声功率

Error_Count_1T1R=zeros(1,length(SNR_dB));  %码元错误个数（1发1收）
Error_Count_2T1R=zeros(1,length(SNR_dB));  %码元错误个数（2发1收）
Error_Count_1T2R=zeros(1,length(SNR_dB));  %码元错误个数（1发2收）

SER_1T1R=zeros(1,length(SNR_dB));          %仿真误码率（1发1收）
SER_2T1R=zeros(1,length(SNR_dB));          %仿真误码率（2发1收）
SER_1T2R=zeros(1,length(SNR_dB));          %仿真误码率（1发2收）

Demod_1T1R=zeros(1,N);                     %解调信号（1发1收）
Demod_2T1R=zeros(2,N);                     %解调信号（2发1收）
Demod_1T2R=zeros(1,N);                     %解调信号（1发2收）

yn_1T1R=zeros(1,N);                        %接收信号（1发1收）
yn_2T1R=zeros(2,N);                        %接收信号（2发1收）
yn_1T2R=zeros(1,N);                        %接收信号（1发2收）

h_1T1R=(randn(1,N) +1i*randn(1,N));        %复高斯瑞利信道,均值为0+0i，方差为2（1发1收）
h_2T1R=(randn(2,N) +1i*randn(2,N));        %复高斯瑞利信道,均值为0+0i，方差为2（2发1收）
h_1T2R=(randn(2,N) +1i*randn(2,N));        %复高斯瑞利信道,均值为0+0i，方差为2（1发2收）

h_2T1R=h_2T1R./sqrt(2);                    %信道能量归一化，方差为1（2发1收） 
h_1T2R=h_1T2R./sqrt(2);                    %信道能量归一化，方差为1（1发2收） 
h_1T1R=h_1T1R./sqrt(2);                    %信道能量归一化，方差为1（1发1收） 

%% 1【1发1收】
sn_1T1R=round(rand(1,N));              % 原信号： 0,1序列sn
sn1_1T1R=round((sn_1T1R-1/2)*2);       % 原信号：-1,1序列sn1
for i=1:length(SNR_dB)
    noise_1T1R = sqrt(Pn(i)/2)*(randn(1,N)+1i*randn(1,N)); %高斯白噪声
%     yn_1T1R=sn1_1T1R+noise_1T1R./h_1T1R; 
    r_1T1R=h_1T1R.*sn1_1T1R+noise_1T1R; %接收
    yn_1T1R=r_1T1R./h_1T1R;  %判决变量
    %%% 解码
    for j=1:N
         if (yn_1T1R(j)>0)
             Demod_1T1R(j)=1;         %接收信号大于0，则判1
         else
             Demod_1T1R(j)=0;         %接收信号小于0，则判0
         end
    end
     %%% 计算错误比特个数
     for j=1:N
         if(Demod_1T1R(j)~=sn_1T1R(j))
              Error_Count_1T1R(i)=Error_Count_1T1R(i)+1;
         end
     end
     %%% 计算误码率
     SER_1T1R(i)=Error_Count_1T1R(i)/N;           % 仿真误码率
end
%% 2【2发1收 空时码】
sn_2T1R=round(rand(2,N));                           % 原信号： 0,1序列sn
sn1_2T1R=(1/sqrt(2)).*round((sn_2T1R-1/2)*2);       % 原信号

for i=1:length(SNR_dB)
    %%% 生成信号
    noise_1= sqrt(Pn(i)/2)*(randn(1,N)+1i*randn(1,N)); %第一个信道上的噪声
    noise_2= sqrt(Pn(i)/2)*(randn(1,N)+1i*randn(1,N)); %第二个信道上的噪声
%     noise_2T1R(1,:)= (conj(h_2T1R(1,:)).*noise_1+h_2T1R(2,:).*conj(noise_2))./(sum(abs(h_2T1R).^2)); %(h1*n1+h2n2*)/(||h||)             
%     noise_2T1R(2,:)= (conj(h_2T1R(2,:)).*noise_1-h_2T1R(1,:).*conj(noise_2))./(sum(abs(h_2T1R).^2)); %(h2*n1-h1n2*)/(||h||)
%     yn_2T1R=sn1_2T1R+noise_2T1R; 
    r1_2T1R=h_2T1R(1,:).*sn1_2T1R(1,:)+h_2T1R(2,:).*sn1_2T1R(2,:)+noise_1; %t1时刻接收
    r2_2T1R=-h_2T1R(1,:).*conj(sn1_2T1R(2,:))+h_2T1R(2,:).*conj(sn1_2T1R(1,:))+noise_2; %t2时刻接收
    yn_2T1R(1,:)=( conj(h_2T1R(1,:)).*r1_2T1R+h_2T1R(2,:).*conj(r2_2T1R) )./(sum(abs(h_2T1R).^2)); %判决变量
    yn_2T1R(2,:)=( conj(h_2T1R(2,:)).*r1_2T1R-h_2T1R(1,:).*conj(r2_2T1R) )./(sum(abs(h_2T1R).^2));
    %%% 解码
    for j=1:2
        for k=1:N
           if (yn_2T1R(j,k)>0)
             Demod_2T1R(j,k)=1;       %接收信号大于0，则判1
           else
             Demod_2T1R(j,k)=0;       %接收信号小于0，则判0
           end
        end
    end
    %%% 计算错误比特个数
    for j=1:2
        for k=1:N
          if(Demod_2T1R(j,k)~=sn_2T1R(j,k))
              Error_Count_2T1R(i)=Error_Count_2T1R(i)+1;
          end
        end
    end
    %%% 计算误码率
    SER_2T1R(i)=Error_Count_2T1R(i)/N/2;           % 仿真误码率
end

%% 3【1发2收】
sn_1T2R=round(rand(1,N));                                       % 原信号： 0,1序列sn
sn1_1T2R=[round((sn_1T2R-1/2)*2);round((sn_1T2R-1/2)*2)];       % 原信号：-1,1序列sn1
for i=1:length(SNR_dB)
    %%% 生成信号
    noise_1T2R=sqrt(Pn(i)/2)*(randn(2,N)+1i*randn(2,N));
    r_1T2R=h_1T2R.*sn1_1T2R+noise_1T2R;  %接收信号
%     yn_1T2R=(conj(h_1T2R(1,:)).*(h_1T2R(1,:).*sn1_1T2R(1,:)+noise_1T2R(1,:))+... %y=[h1*(h1x1+n1)+h2*(h2x2+n2)]/(h1^2+h2^2)
%              conj(h_1T2R(2,:)).*(h_1T2R(2,:).*sn1_1T2R(2,:)+noise_1T2R(2,:)))./...
%              (abs(h_1T2R(1,:)).^2+abs(h_1T2R(2,:)).^2);
    yn_1T2R=(conj(h_1T2R(1,:)).*r_1T2R(1,:)+conj(h_1T2R(2,:)).*r_1T2R(2,:))./(sum(abs(h_1T2R).^2));
    %%% 解码
    for j=1:N
        if (yn_1T2R(j)>0)
             Demod_1T2R(j)=1;       %接收信号大于0，则判1
        else
             Demod_1T2R(j)=0;       %接收信号小于0，则判0
        end
    end
    %%% 计算错误比特个数
    for j=1:N
       if(Demod_1T2R(j)~=sn_1T2R(j))
           Error_Count_1T2R(i)=Error_Count_1T2R(i)+1;
       end
    end
    %%% 计算误码率
    SER_1T2R(i)=Error_Count_1T2R(i)/N;           % 仿真误码率
end

%% [4、显示结果]
semilogy(SNR_dB,SER_1T1R,'go-',SNR_dB,SER_2T1R,'r*-',SNR_dB,SER_1T2R,'k*-');     
axis([0 45 10^-5 10^0]);     
legend("1发1收","2发1收","1发2收");
xlabel('信噪比SNR/dB');  %横轴坐标
ylabel('误码率BER');     %纵轴坐标
title("BPSK调制信号在复高斯瑞利信道下Alamouti的误码率分析");
grid on;                    %显示表格线
%     yn_1T2R=( r_1T2R*conj(h_1T2R'))./(sum(abs(h_1T2R).^2)); %判决变量./(sum(abs(h_1T2R).^2))