function BER = Alamouti_BPSK( num_trans, num_rx, SNR)
%計算BPSK調制方式下Alamouti編碼的誤比特率
%BER：誤比特率；num_trans：發送的比特組數；
%num_rx：接收天線的個數；SNR：給定的信噪比
%

if (num_trans < 1) || (num_rx < 1) || (SNR < 0)
    disp('參數錯誤');
    return;
end;

num_biterr_total = 0;   %總的誤比特個數
for loop = 1 : num_trans
    
    %產生信源
    %S：2×2，sij：第i根發射天線在第j個slot發送的信號
    S = zeros(2, 2);
    seed = sign(randn(1, 2));
    S(1, 1) = seed(1);
    S(2, 1) = seed(2);
    S(1, 2) = -conj(seed(2));
    S(2, 2) = conj(seed(1));
    
    %產生信道衰落系數
    %H:2×num_rx，hij:第i根發射天線和第j根接收天線之間的衰落系數
    H = zeros(2, num_rx);
    A = randn(2, num_rx);
    B = randn(2, num_rx);
    H = complex(A, B);
    H = H / sqrt(2);    %信道矩陣歸一化
    
    %產生接收端AWGN
    %N:2×num_rx，nij：第i個slot內第j根接收天線處的噪聲
    sigma = 2 * 10 ^ (-SNR / 10);   %噪聲的功率；因為有兩根發射天線，每根發射天線的發射功率均為1，所以要乘以2
    N = zeros(2, num_rx);
    C = randn(2, num_rx);
    D = randn(2, num_rx);
    N = complex(C, D);
    N = N / sqrt(2);
    N = sqrt(sigma) * N;    %噪聲矩陣歸一化，使其中每一個元素的功率＝sigma
    
    %接收信號
    %產生接收矩陣：R:2×num_rx,rij:第i個slot內第j根接收天線接收到的信號
    R = zeros(2, num_rx);
    R = S' * H + N;  %生成接收信號
    
    %構造判決統計
    X_kepa = zeros(2, 1);   %最終判決統計結果
    temp1 = zeros(1, num_rx);
    temp2 = zeros(1, num_rx);
    temp3 = zeros(1, num_rx);
    temp4 = zeros(1, num_rx);
    
    temp1 = R(1, :) .* conj(H(1, :));
    temp2 = conj(R(2, :)) .* H(2, :);
    temp3 = R(1, :) .* conj(H(2, :));
    temp4 = conj(R(2, :)) .* H(1, :);
    
    deposit1 = 0;
    deposit2 = 0;
    deposit3 = 0;
    deposit4 = 0;
    
    for subloop = 1 : num_rx
        deposit1 = deposit1 + temp1(subloop);
        deposit2 = deposit2 + temp2(subloop);
        deposit3 = deposit3 + temp3(subloop);
        deposit4 = deposit4 + temp4(subloop);
    end
    
    X_kepa(1, 1) = deposit1 + deposit2;
    X_kepa(2, 1) = deposit3 - deposit4;
    
    %進行最大似然譯碼
    dist = zeros(2, 2); %距離矩陣
    result = zeros(1, 2); %譯碼結果
    
    dist(1, 1) = Eucl_Dist(X_kepa(1, 1), -1);
    dist(1, 2) = Eucl_Dist(X_kepa(1, 1), 1);
    dist(2, 1) = Eucl_Dist(X_kepa(2, 1), -1);
    dist(2, 2) = Eucl_Dist(X_kepa(2, 1), 1);
    
    if dist(1, 1) < dist(1, 2)
        result(1) = -1;
    else
        result(1) = 1;
    end;
    if dist(2, 1) < dist(2, 2)
        result(2) = -1;
    else
        result(2) = 1;
    end;
    
    %統計誤比特個數
    num_biterr_temp = 0;
    if result(1) ~= S(1, 1)
        num_biterr_temp = num_biterr_temp + 1;
    end;
    if result(2) ~= S(2, 1)
        num_biterr_temp = num_biterr_temp + 1;
    end;
    num_biterr_total = num_biterr_total + num_biterr_temp;
end
%計算誤比特率
BER = num_biterr_total / (num_trans * 2);