% 對BCH碼進行了仿真，可變調制方式BPSK,4QAM,16QAM等
% 每楨長度等于K的整數倍
clear 
n_frame=20;
N=15;
K=11;
M=1;    %1--BPSK    2--QPSK or QAM    3---8PSK   4---16PSK or 16QAM
          % psk        psk     qask        psk         psk     qask
Es=1;   % 數字映射后的符號能量 psk--1
        %                      16QAM--10
L_frame=10000;
start=-16;
step=2;
finish=20;

[pg, pm, cs, h, t] = bchpoly(N, K);
r=K/N; % 編碼效率
Eb=Es/M; % 每比特能量
plot_pe=[];
Q=1;
axis_EbN0=start:step:finish;
%******************************************************************8

for EbN0=start:step:finish
    Liner_EbN0=10^(EbN0/10);
    pe_number=0;
    variance=0.5*(Eb/Liner_EbN0)/r;  % 噪聲方差, 來源于a^2/(2*delta^2)=Liner_EbN0,a^2信息幅值的平方即比特能量, 2*delta^2=N0.    N0/r得到編碼后的噪聲功率譜密度
    for i=1:1:n_frame
        
        x_msg=randint(L_frame/K,K,2);     %真正的信息
       
        x_code_msg=bchenco(x_msg,N,K);    % 編碼模塊
        
        X_code_msg=reshape(x_code_msg',1,[]); % 碼字整形，矩陣變成行向量
       
        symbol=bi2symbol(X_code_msg,M); % 二進制到十進制轉換
        
        modulation_sig=modmap(symbol,1,1,'qask',2^M); % 數字映射
        
        x_ch=g_ch(modulation_sig,variance);  % 經過高斯信道，噪聲方差var
        
        de_symbol=demodmap(x_ch,1,1,'qask',2^M); %逆映射
        
        rec=symbol2bit(de_symbol,M); % 十進制轉換二進制
        
        
        
        out=reshape(rec,N,L_frame/K)'; % 碼字整形，行向量變成矩陣
        in=out;
        
        x_dec_msg=bchdeco(in,K,t);% 解碼模塊
        
        pe_number=pe_number+sum( sum(x_msg~=x_dec_msg) );
       
        current_time=fix(clock);
       fprintf('i am working %g      %g年  %g月  %g日  %g時  %g分  %g秒\n\n',Q,current_time(1),current_time(2),current_time(3),current_time(4),current_time(5),current_time(6))
        Q=Q+1;
        fprintf('\n\n')
        
    end
    pe=pe_number/(L_frame*n_frame);
    plot_pe=[plot_pe,pe];
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
semilogy(axis_EbN0,plot_pe,'b*-')

xlabel('Eb/N0 信號功率/噪聲功率(db)')
ylabel('BER')
title('誤碼率')
grid on