ml_and_map.m

最大似然（ML）准则和最大后验概率（MAP）准则Matlab仿真

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准则信噪比与无码率的曲线