%%16PSK的調制與解調,programmed by阿哲
clear all
clc
N = 2*10^3;
M = 16;
K = 4;
Err_bits = 0;

%%16psk modu
Tx_values = rand(1,N)-0.5;
for i=1:length(Tx_values)
    if Tx_values(i)>=0
        Tx_bits(i)=1;
    else
        Tx_bits(i)=0;
    end
end

for i=1:4:length(Tx_bits)
    bin_data = Tx_bits(i:i+3);
    dec_data = (8*Tx_bits(i))+(4*Tx_bits(i+1))+(2*Tx_bits(i+2))+Tx_bits(i+3);
    Mpsk_modu((i-1)/4+1) = exp(j*dec_data*2*pi/M);
    dec_data = 0;
end


%%IQ modu
Mpsk_modu_real = real(Mpsk_modu);
Mpsk_modu_imag = imag(Mpsk_modu);

scatterplot(Mpsk_modu)
title('16PSK信號的星座圖(發送端)')

figure
subplot(211)
plot(Mpsk_modu_real)
xlabel('同相路信號')
ylabel('幅度')
title('發送端的16PSK信號')
subplot(212)
plot(Mpsk_modu_imag)
xlabel('正交路信號')
ylabel('幅度')
title('發送端的16PSK信號')

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%Channel
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Mpsk_modu = awgn(complex(Mpsk_modu_real,Mpsk_modu_imag),25);

scatterplot(Mpsk_modu)
title('16PSK信號的星座圖（接收端）')

%%16psk demodu
for j = 1:length(Mpsk_modu)
    ang_Mpsk(j) = angle(Mpsk_modu(j));
    if ang_Mpsk(j)>=pi/16&&ang_Mpsk(j)<pi*3/16
        Demapper_output((j-1)*K+1) = 0;
        Demapper_output((j-1)*K+2) = 0;
        Demapper_output((j-1)*K+3) = 0;
        Demapper_output((j-1)*K+4) = 1;
    elseif ang_Mpsk(j)>=pi*3/16&&ang_Mpsk(j)<pi*5/16
        Demapper_output((j-1)*K+1) = 0;
        Demapper_output((j-1)*K+2) = 0;
        Demapper_output((j-1)*K+3) = 1;
        Demapper_output((j-1)*K+4) = 0;
    elseif ang_Mpsk(j)>=pi*5/16&&ang_Mpsk(j)<pi*7/16
        Demapper_output((j-1)*K+1) = 0;
        Demapper_output((j-1)*K+2) = 0;
        Demapper_output((j-1)*K+3) = 1;
        Demapper_output((j-1)*K+4) = 1;
    elseif ang_Mpsk(j)>=pi*7/16&&ang_Mpsk(j)<pi*9/16
        Demapper_output((j-1)*K+1) = 0;
        Demapper_output((j-1)*K+2) = 1;
        Demapper_output((j-1)*K+3) = 0;
        Demapper_output((j-1)*K+4) = 0;
    elseif ang_Mpsk(j)>=pi*9/16&&ang_Mpsk(j)<pi*11/16
        Demapper_output((j-1)*K+1) = 0;
        Demapper_output((j-1)*K+2) = 1;
        Demapper_output((j-1)*K+3) = 0;
        Demapper_output((j-1)*K+4) = 1;
    elseif ang_Mpsk(j)>=pi*11/16&&ang_Mpsk(j)<pi*13/16
        Demapper_output((j-1)*K+1) = 0;
        Demapper_output((j-1)*K+2) = 1;
        Demapper_output((j-1)*K+3) = 1;
        Demapper_output((j-1)*K+4) = 0;
    elseif ang_Mpsk(j)>=pi*13/16&&ang_Mpsk(j)<pi*15/16
        Demapper_output((j-1)*K+1) = 0;
        Demapper_output((j-1)*K+2) = 1;
        Demapper_output((j-1)*K+3) = 1;
        Demapper_output((j-1)*K+4) = 1;
    elseif ang_Mpsk(j)>=pi*15/16||ang_Mpsk(j)<-pi*15/16
        Demapper_output((j-1)*K+1) = 1;
        Demapper_output((j-1)*K+2) = 0;
        Demapper_output((j-1)*K+3) = 0;
        Demapper_output((j-1)*K+4) = 0;
    elseif ang_Mpsk(j)>=-pi*15/16&&ang_Mpsk(j)<-pi*13/16
        Demapper_output((j-1)*K+1) = 1;
        Demapper_output((j-1)*K+2) = 0;
        Demapper_output((j-1)*K+3) = 0;
        Demapper_output((j-1)*K+4) = 1;
    elseif ang_Mpsk(j)>=-pi*13/16&&ang_Mpsk(j)<-pi*11/16
        Demapper_output((j-1)*K+1) = 1;
        Demapper_output((j-1)*K+2) = 0;
        Demapper_output((j-1)*K+3) = 1;
        Demapper_output((j-1)*K+4) = 0;
    elseif ang_Mpsk(j)>=-pi*11/16&&ang_Mpsk(j)<-pi*9/16
        Demapper_output((j-1)*K+1) = 1;
        Demapper_output((j-1)*K+2) = 0;
        Demapper_output((j-1)*K+3) = 1;
        Demapper_output((j-1)*K+4) = 1;
    elseif ang_Mpsk(j)>=-pi*9/16&&ang_Mpsk(j)<-pi*7/16
        Demapper_output((j-1)*K+1) = 1;
        Demapper_output((j-1)*K+2) = 1;
        Demapper_output((j-1)*K+3) = 0;
        Demapper_output((j-1)*K+4) = 0;
    elseif ang_Mpsk(j)>=-pi*7/16&&ang_Mpsk(j)<-pi*5/16
        Demapper_output((j-1)*K+1) = 1;
        Demapper_output((j-1)*K+2) = 1;
        Demapper_output((j-1)*K+3) = 0;
        Demapper_output((j-1)*K+4) = 1;
    elseif ang_Mpsk(j)>=-pi*5/16&&ang_Mpsk(j)<-pi*3/16
        Demapper_output((j-1)*K+1) = 1;
        Demapper_output((j-1)*K+2) = 1;
        Demapper_output((j-1)*K+3) = 1;
        Demapper_output((j-1)*K+4) = 0;
    elseif ang_Mpsk(j)>=-pi*3/16&&ang_Mpsk(j)<-pi*1/16
        Demapper_output((j-1)*K+1) = 1;
        Demapper_output((j-1)*K+2) = 1;
        Demapper_output((j-1)*K+3) = 1;
        Demapper_output((j-1)*K+4) = 1;
    elseif ang_Mpsk(j)>=-pi*1/16&&ang_Mpsk(j)<pi*1/16
        Demapper_output((j-1)*K+1) = 0;
        Demapper_output((j-1)*K+2) = 0;
        Demapper_output((j-1)*K+3) = 0;
        Demapper_output((j-1)*K+4) = 0;
    end
end

Rv_bits = Demapper_output;

%%Compare BER
for l = 1:length(Rv_bits)
    if Rv_bits(l)~=Tx_bits(l)
        Err_bits = Err_bits+1;
    end
end