【MATLAB】各種変調方式のAWGN通信でのBER特性

以下のコードを実行すれば各種変調方式(BPSK,QPSK,16QAM,64QAM)におけるBER特性を求められるはずです。

clear all

Nt = 2; %送信アンテナ数
L = 6 * 10^6; %送信ビット数
EbNo = 0:16;

RX_bpsk_symbol = zeros(Nt,L);
RX_bpsk_bit = zeros(Nt,L);

BER = zeros(4,length(EbNo));

for count_s = 1:length(EbNo)
    TX_bit = randi([0 1],Nt,L);

    TX_bpsk = ones(Nt,L);
    TX_bpsk(TX_bit == 0) = -1;

    TX_qpsk = bit_2_QPSK(TX_bpsk);
    TX_16qam = bit_2_16QAM(TX_bpsk);
    TX_64qam = bit_2_64QAM(TX_bpsk);

    RX_bpsk = AWGNchannel(TX_bpsk,EbNo(count_s),1);
    RX_qpsk = AWGNchannel(TX_qpsk,EbNo(count_s),2);
    RX_16qam = AWGNchannel(TX_16qam,EbNo(count_s),4);
    RX_64qam = AWGNchannel(TX_64qam,EbNo(count_s),6);

    RX_bpsk_symbol(real(RX_bpsk) < 0) = -1;
    RX_bpsk_symbol(real(RX_bpsk) > 0) = 1;

    RX_qpsk_symbol = QPSK_2_symbol(RX_qpsk);
    RX_16qam_symbol = QAM16_2_symbol(RX_16qam);
    RX_64qam_symbol = QAM64_2_symbol(RX_64qam);

    RX_bpsk_bit(RX_bpsk_symbol == 1) = 1;
    RX_bpsk_bit(RX_bpsk_symbol == -1) = 0;

    RX_qpsk_bit = symbol_2_bit_QPSK(RX_qpsk_symbol);
    RX_16qam_bit = symbol_2_bit_16QAM(RX_16qam_symbol);
    RX_64qam_bit = symbol_2_bit_64QAM(RX_64qam_symbol);

    RX_qpsk_bit(RX_qpsk_bit == 1) = 1;
    RX_qpsk_bit(RX_qpsk_bit == -1) = 0;
    RX_16qam_bit(RX_16qam_bit == 1) = 1;
    RX_16qam_bit(RX_16qam_bit == -1) = 0;
    RX_64qam_bit(RX_64qam_bit == 1) = 1;
    RX_64qam_bit(RX_64qam_bit == -1) = 0;

    BER(1,count_s) = sum(abs(TX_bit(:) - RX_bpsk_bit(:)))/(Nt*L);
    BER(2,count_s) = sum(abs(TX_bit(:) - RX_qpsk_bit(:)))/(Nt*L);
    BER(3,count_s) = sum(abs(TX_bit(:) - RX_16qam_bit(:)))/(Nt*L);
    BER(4,count_s) = sum(abs(TX_bit(:) - RX_64qam_bit(:)))/(Nt*L);
    disp(['Eb/No = ' num2str(EbNo(count_s)) ' finished !'])
end

QPSK関連の関数の中身はこちら：https://qiita.com/Seiji_Tanaka/items/f8958c6d2f4b4780de5b
16QAM関連の中身はこちら：https://qiita.com/Seiji_Tanaka/items/add956f5a692bc18cf65
64QAM関連の中身はこちら：https://qiita.com/Seiji_Tanaka/items/2a3c822bdfea6dcf5412
AWGNchannelの中身はこちら：https://qiita.com/Seiji_Tanaka/items/f033e56b54e596637bfa

グラフ出力用のコードも貼っておきます。

figure

semilogy(EbNo,BER(2,:),'-^','Color',[252,200,0]./255,'LineWidth',3,'MarkerSize',12,'MarkerFaceColor',[1,1,1])
 hold on
semilogy(EbNo,BER(3,:),'-s','Color',[0,145,197]./255,'LineWidth',3,'MarkerSize',12,'MarkerFaceColor',[1,1,1])
hold on
semilogy(EbNo,BER(4,:),'-v','Color',[40,182,89]./255,'LineWidth',3,'MarkerSize',12,'MarkerFaceColor',[1,1,1])
grid on

set(gca,'FontSize',20);

ax = gca;
ax.LineWidth = 2;

set( gca, 'FontName','Times','FontSize',20 ); 

ylim([10^(-6) 10^(-0)])
xlim([EbNo(1) EbNo(end)])

xlabel('{\sl E}_{\rm b}/{\sl N}_{\rm 0} [dB]','FontSize',20,'Interpreter','tex')
ylabel('BER','FontSize',20)
lgd = legend({'QPSK','16QAM','64QAM'},'Location', 'best');
lgd.FontSize = 12;

こんな感じで結果がグラフが出力されます。