function [llr]=deqam16(x,h,SNR)
%本程序是計算經過16QAM調制信號的對數釋然比信息，為軟維特比譯碼提供信息源
%這里是根據輸入信息點在坐標中的位置情況來判斷對應二進制信息為的似然比的大小
%這里應該是對應四位二進制比特的釋然比
%根據二進制對應比特數位置，分別求離輸入信息點最近的0和1對應星座的似然比即位該比特的似然比
%具體理論公式見第一階段文檔附件中的均衡和解調方案
%輸入參數：x：  經過16QAM調制的信號，長度為L
%          h：  信號經過的信號估計序列
%          SNR：經過的高斯信道的信噪比
%輸出參數：llr：信號對應二進制序列的對數釋然比信息
%                    8    7  |   3    4
%                    6    5  |   1    2
%             －－－－－－－－－－－－－－－        
%                   14    13 |   9    10
%                   16    15 |   11   12
%
%王東洋
%2005－10－19
%星座圖表
constel_diagram=[sqrt(2)/2+j*sqrt(2)/2, 1.5*sqrt(2)+j*sqrt(2)/2, sqrt(2)/2+j*1.5*sqrt(2), 1.5*sqrt(2)+j*1.5*sqrt(2),...%第一象限
                -sqrt(2)/2+j*sqrt(2)/2,-1.5*sqrt(2)+j*sqrt(2)/2,-sqrt(2)/2+j*1.5*sqrt(2),-1.5*sqrt(2)+j*1.5*sqrt(2),...%第二象限
                 sqrt(2)/2-j*sqrt(2)/2, 1.5*sqrt(2)-j*sqrt(2)/2, sqrt(2)/2-j*1.5*sqrt(2), 1.5*sqrt(2)-j*1.5*sqrt(2),...%第四象限
                -sqrt(2)/2-j*sqrt(2)/2,-1.5*sqrt(2)-j*sqrt(2)/2,-sqrt(2)/2-j*1.5*sqrt(2),-1.5*sqrt(2)-j*1.5*sqrt(2)]/sqrt(5);%第三象限
%SNR為信道信噪比，信號功率為1,轉化為線性值
SNR_linr=10^(SNR/10);
%得到信道估計值的平方
h_square=abs(h).^2;
%len為輸入信號的長度
len=length(x);
%存儲似然比信息
llr=zeros(1,4*len);
%計算信號到各個星座點的映射距離
temp=[abs(x-constel_diagram(1)),abs(x-constel_diagram(2)),abs(x-constel_diagram(3)),abs(x-constel_diagram(4)),abs(x-constel_diagram(5)),abs(x-constel_diagram(6)),...
    abs(x-constel_diagram(7)),abs(x-constel_diagram(8)),abs(x-constel_diagram(9)),abs(x-constel_diagram(10)),abs(x-constel_diagram(11)),abs(x-constel_diagram(12)),...
    abs(x-constel_diagram(13)),abs(x-constel_diagram(14)),abs(x-constel_diagram(15)),abs(x-constel_diagram(16))].^2;
%下面是計算各個信息比特的似然比信息llr
for m=1:len
    temp2=temp(m:len:16*len);
    [y,z]=sort(temp2);
    
    switch z(1)
        case 1
        llr(4*(m-1)+1)=h_square(m)*(temp2(5)-temp2(1));
        llr(4*(m-1)+2)=h_square(m)*(temp2(9)-temp2(1));
        llr(4*(m-1)+3)=h_square(m)*(temp2(2)-temp2(1));
        llr(4*(m-1)+4)=h_square(m)*(temp2(3)-temp2(1));
        case 2
         llr(4*(m-1)+1)=h_square(m)*(temp2(5)-temp2(2));
         llr(4*(m-1)+2)=h_square(m)*(temp2(10)-temp2(2));
         llr(4*(m-1)+3)=h_square(m)*(temp2(2)-temp2(1));
         llr(4*(m-1)+4)=h_square(m)*(temp2(4)-temp2(2));
        case 3
         llr(4*(m-1)+1)=h_square(m)*(temp2(7)-temp2(3));
         llr(4*(m-1)+2)=h_square(m)*(temp2(9)-temp2(3));
         llr(4*(m-1)+3)=h_square(m)*(temp2(4)-temp2(3));
         llr(4*(m-1)+4)=h_square(m)*(temp2(3)-temp2(1));
        case 4
         llr(4*(m-1)+1)=h_square(m)*(temp2(7)-temp2(4));
         llr(4*(m-1)+2)=h_square(m)*(temp2(10)-temp2(4));
         llr(4*(m-1)+3)=h_square(m)*(temp2(4)-temp2(3));
         llr(4*(m-1)+4)=h_square(m)*(temp2(4)-temp2(2));
        case 5
         llr(4*(m-1)+1)=h_square(m)*(temp2(5)-temp2(1));
         llr(4*(m-1)+2)=h_square(m)*(temp2(13)-temp2(5));
         llr(4*(m-1)+3)=h_square(m)*(temp2(6)-temp2(5));
         llr(4*(m-1)+4)=h_square(m)*(temp2(7)-temp2(5));
        case 6
         llr(4*(m-1)+1)=h_square(m)*(temp2(6)-temp2(1));
         llr(4*(m-1)+2)=h_square(m)*(temp2(14)-temp2(6));
         llr(4*(m-1)+3)=h_square(m)*(temp2(6)-temp2(5));
         llr(4*(m-1)+4)=h_square(m)*(temp2(8)-temp2(6));
        case 7
         llr(4*(m-1)+1)=h_square(m)*(temp2(7)-temp2(3));
         llr(4*(m-1)+2)=h_square(m)*(temp2(13)-temp2(7));
         llr(4*(m-1)+3)=h_square(m)*(temp2(8)-temp2(7));
         llr(4*(m-1)+4)=h_square(m)*(temp2(7)-temp2(5));
        case 8
         llr(4*(m-1)+1)=h_square(m)*(temp2(8)-temp2(3));
         llr(4*(m-1)+2)=h_square(m)*(temp2(14)-temp2(6));
         llr(4*(m-1)+3)=h_square(m)*(temp2(8)-temp2(7));
         llr(4*(m-1)+4)=h_square(m)*(temp2(8)-temp2(6));
        case 9
         llr(4*(m-1)+1)=h_square(m)*(temp2(13)-temp2(9));
         llr(4*(m-1)+2)=h_square(m)*(temp2(9)-temp2(1));
         llr(4*(m-1)+3)=h_square(m)*(temp2(10)-temp2(9));
         llr(4*(m-1)+4)=h_square(m)*(temp2(11)-temp2(9));
        case 10
         llr(4*(m-1)+1)=h_square(m)*(temp2(13)-temp2(10));
         llr(4*(m-1)+2)=h_square(m)*(temp2(10)-temp2(2));
         llr(4*(m-1)+3)=h_square(m)*(temp2(10)-temp2(9));
         llr(4*(m-1)+4)=h_square(m)*(temp2(12)-temp2(10));
        case 11
         llr(4*(m-1)+1)=h_square(m)*(temp2(15)-temp2(11));
         llr(4*(m-1)+2)=h_square(m)*(temp2(11)-temp2(1));
         llr(4*(m-1)+3)=h_square(m)*(temp2(12)-temp2(11));
         llr(4*(m-1)+4)=h_square(m)*(temp2(11)-temp2(9));
        case 12
         llr(4*(m-1)+1)=h_square(m)*(temp2(15)-temp2(12));
         llr(4*(m-1)+2)=h_square(m)*(temp2(12)-temp2(2));
         llr(4*(m-1)+3)=h_square(m)*(temp2(12)-temp2(11));
         llr(4*(m-1)+4)=h_square(m)*(temp2(12)-temp2(10));
        case 13
         llr(4*(m-1)+1)=h_square(m)*(temp2(13)-temp2(9));
         llr(4*(m-1)+2)=h_square(m)*(temp2(13)-temp2(5));
         llr(4*(m-1)+3)=h_square(m)*(temp2(14)-temp2(13));
         llr(4*(m-1)+4)=h_square(m)*(temp2(15)-temp2(13));
        case 14
         llr(4*(m-1)+1)=h_square(m)*(temp2(14)-temp2(9));
         llr(4*(m-1)+2)=h_square(m)*(temp2(14)-temp2(6));
         llr(4*(m-1)+3)=h_square(m)*(temp2(14)-temp2(13));
         llr(4*(m-1)+4)=h_square(m)*(temp2(16)-temp2(14));
        case 15
         llr(4*(m-1)+1)=h_square(m)*(temp2(15)-temp2(11));
         llr(4*(m-1)+2)=h_square(m)*(temp2(15)-temp2(5));
         llr(4*(m-1)+3)=h_square(m)*(temp2(16)-temp2(15));
         llr(4*(m-1)+4)=h_square(m)*(temp2(15)-temp2(13));
       otherwise 
         llr(4*(m-1)+1)=h_square(m)*(temp2(16)-temp2(11));
         llr(4*(m-1)+2)=h_square(m)*(temp2(16)-temp2(6));
         llr(4*(m-1)+3)=h_square(m)*(temp2(16)-temp2(15));
         llr(4*(m-1)+4)=h_square(m)*(temp2(16)-temp2(14));
     end
 end
 llr=llr/SNR_linr;