clear
clc
%輸入參數變量的值
p = 0.1;    %序列{+1,-1}中+1出現的概率為p，－1出現的概率為1-p
signal_power = 1;   %信號功率為1
noise_power = 0.0;   %噪聲功率為0
noise_power_times = 50; %信噪比的次數
test_times = 100;  %重復實驗的次數
length = 1000;  %序列的長度
SNR_dB = zeros(noise_power_times,1);  %保留每次實驗的SNR
BER_mean_MAP = zeros(noise_power_times,1);  %保留每次MAP實驗BER的均值
BER_mean_ML = zeros(noise_power_times,1);  %保留每次ML實驗BER的均值
BER_all_MAP = zeros(test_times,1);  %保留每次MAP實驗的BER
BER_all_ML = zeros(test_times,1);  %保留每次ML實驗的BER
%運算過程
for i = 1: noise_power_times    %信噪比的個數
    noise_power = noise_power+i/10; %噪聲功率N
    noise_mean = 0; %高斯白噪聲的均值
    noise_var = noise_power*0.5;    %高斯白噪聲的方差    
    for j = 1:test_times    %進行重復實驗
        %信號序列
        signal = randsrc(1,length,[-1 1;1-p p]);   %randsrc()產生-1和1的貝努利隨機序列，概率均為 0.5
        %噪聲序列
        noise = normrnd(noise_mean,noise_var,[1,length]);   %normrnd()產生N(noise_mean,noise_var)的高斯白噪聲
        %加噪聲的信號序列
        mixed_signal = signal + noise;
        judged_signal_MAP = zeros(size(signal));   %zeros()初始化序列的值為0，size()獲得矩陣的行數和列數
        judged_signal_ML = zeros(size(signal));   %zeros()初始化序列的值為0，size()獲得矩陣的行數和列數
        for k = 1:length
            %最佳MAP判決的信號序列
            if mixed_signal(k) >= 0.25*noise_power*log((1-p)/p)   %判決門限
                 judged_signal_MAP(k) = 1;
            else                      
                judged_signal_MAP(k) = -1;
            end   
            %最佳ML判決的信號序列
            if mixed_signal(k) >= 0   %判決門限
                 judged_signal_ML(k) = 1;
            else                      
                judged_signal_ML(k) = -1;
            end   
        end 
        %計算兩個序列的差別，判決的信號和原來的信號相同的位置為 0
        difference_MAP = signal - judged_signal_MAP;
        difference_ML = signal - judged_signal_ML;
        %計算錯誤bit的個數
        bit_error_MAP = nnz(difference_MAP);   %nnz()計算非0值的個數
        bit_error_ML = nnz(difference_ML);
        %計算誤比特率
        BER_all_MAP(j) = bit_error_MAP / length;
        BER_all_ML(j) = bit_error_ML / length;
    end   %對應for i = 1:times

    %計算誤比特率的數學期望
    BER_mean_MAP(i) = mean(BER_all_MAP);   %mean()計算平均值
    BER_mean_ML(i) = mean(BER_all_ML);   %mean()計算平均值
    SNR_dB(i) = 10*LOG10(signal_power/noise_power);
end
%繪出圖形
figure(1);   %figure()新建一個繪圖窗口
plot(SNR_dB,BER_mean_MAP,'b-o');    %藍線表示最佳MAP準則信噪比與無碼率的曲線
hold on;
plot(SNR_dB,BER_mean_ML,'r-*');     %紅線表示最佳MAP準則信噪比與無碼率的曲線