% 16qam and 64qam ; t4*r4系統; 碼率= 3/4;
% 軟判決譯碼
% 添加定點dsp仿真

% 定點仿真結果： 
% 輸入進行8bits ，16bits輸出，幾乎和16bits量化效果一樣
% 輸入進行8bits ，15bits輸出，幾乎和16bits量化效果一樣
clear
clc
color_flag = 2;
snr_min=0;
snr_max=35;
stepsize = 2.5;
%------------------------仿真參數設定-------------------------------------%
bits = 100; %一個ofdm符號長度
Mod_type = 4; %qam調制 4=16qam;
Max_bit_num = 8;
scale_eff = 3; % 加法器縮放系數
global Max_value Min_value Sub_max_value Sub_min_value
Max_value = 2^(Max_bit_num-1); % Max_value = 2^15
Min_value = -Max_value;
Sub_max_value = 2^(Max_bit_num-1);
Sub_min_value = -Sub_max_value;
Max_add_value = 2^(2*Max_bit_num-1);

Max_div_value = 2^(2*Max_bit_num-3);

mult_flow_num = 0;
add_flow_num = 0;

star_index_length = 16;
Tx_num = 4; %發送天線個數（4）
Rx_num = 4;
Time_length = 4;
pqam_length = bits*Mod_type;
Source_length = (Tx_num-1)*bits*Mod_type; %因為天線3和天線4是一樣的碼子
p_s_length = bits*Mod_type*Time_length;

%-------------------------------------------------------------------------%
j=sqrt(-1);
st2 = sqrt(2);
ch_scale = sqrt(4); % Tx_num = 4
st2 = fix(st2*Max_value)/Max_value;
ch_scale = fix(ch_scale*Max_value)/Max_value;

%------------------產生信號源------------------------%
sig_code = zeros(Tx_num,pqam_length);
sig_code = randint(Tx_num-1,pqam_length);
sig_code(4,:) = sig_code(3,:); 

%-------------------------QAM 映射,STBC編碼--------------%
table16 = [1+j, 3+ j, 1+3j, 3+3j,...
1-j, 1-3j, 3- j, 3-3j,...
-1+j,-1+3j,-3+ j,-3+3j,...
-1-j,-3- j,-1-3j,-3-3j];
table64 = [1+j, 3+ j, 1+3j, 3+3j, 7+1j, 5+ j, 7+3j, 5+3j, 1+7j, 3+7j, 1+5j, 3+5j, 7+7j, 5+7j, 7+5j, 5+5j,...
1-j, 1-3j, 3- j, 3-3j, 1-7j, 1-5j, 3-7j, 3-5j, 7- j, 7-3j, 5- j, 5-3j, 7-7j, 7-5j, 5-7j, 5-5j,...
-1+j,-1+3j,-3+ j,-3+3j,-1+7j,-1+5j,-3+7j,-3+5j,-7+ j,-7+3j,-5+ j,-5+3j,-7+7j,-7+5j,-5+7j,-5+5j,...
-1-j,-3- j,-1-3j,-3-3j,-7- j,-5- j,-7-3j,-5-3j,-1-7j,-3-7j,-1-5j,-3-5j,-7-7j,-5-7j,-7-5j,-5-5j];
table16 = table16/sqrt(10);
table64 = table64/sqrt(42);
table16 = round(table16*Max_value)/Max_value;
table64 = round(table64*Max_value)/Max_value;

qam_code = zeros(Tx_num,bits); %qam映射矩陣的大小(stbc編碼之前)
switch Mod_type
case 4
for i=1:Tx_num
for n=1:bits
table_index = sig_code(i,(n-1)*Mod_type+1)*8+sig_code(i,(n-1)*Mod_type+2)*4+...
sig_code(i,(n-1)*Mod_type+3)*2+sig_code(i,(n-1)*Mod_type+4);
qam_code(i,n) = table16(table_index+1);
end
end
case 6
for i=1:Tx_num
for n=1:bits
table_index = sig_code(i,(n-1)*Mod_type+1)*32+sig_code(i,(n-1)*Mod_type+2)*16+...
sig_code(i,(n-1)*Mod_type+3)*8+sig_code(i,(n-1)*Mod_type+4)*4+...
sig_code(i,(n-1)*Mod_type+5)*2+sig_code(i,(n-1)*Mod_type+6); 
qam_code(i,n) = table64(table_index+1);
end
end 
end
%qam_code = qam_code/sqrt(10);
stbc_code = zeros(Time_length,Tx_num,bits); %stbc 編碼矩陣
stbc_code(1,1,:) = qam_code(1,:);
stbc_code(1,2,:) = qam_code(2,:);
stbc_code(1,3,:) = qam_code(3,:)/st2;
stbc_code(1,4,:) = qam_code(3,:)/st2;

stbc_code(2,1,:) = -conj(qam_code(2,:));
stbc_code(2,2,:) = conj(qam_code(1,:));
stbc_code(2,3,:) = qam_code(3,:)/st2;
stbc_code(2,4,:) = -qam_code(3,:)/st2;

stbc_code(3,1,:) = conj(qam_code(3,:)/st2);
stbc_code(3,2,:) = conj(qam_code(3,:)/st2);
stbc_code(3,3,:) = -real(qam_code(1,:))+j*imag(qam_code(2,:));
stbc_code(3,4,:) = -real(qam_code(2,:))+j*imag(qam_code(1,:));

stbc_code(4,1,:) = conj(qam_code(3,:)/st2);
stbc_code(4,2,:) = -conj(qam_code(3,:)/st2);
stbc_code(4,3,:) = real(qam_code(2,:))+j*imag(qam_code(1,:));
stbc_code(4,4,:) = -(real(qam_code(1,:))+j*imag(qam_code(2,:)));

stbc_code = stbc_code/ch_scale; % 所有發送天線發送信號歸一化


%-----------------------迭代開始 信噪比 由小到大-----------------------------------%
table16 = fix(table16*Max_value);
table64 = fix(table64*Max_value*127/138);
%__________________________________________________________________________________%

iterative_length = 0;
for snr_db_step =snr_min:stepsize:snr_max % 信噪比從snr_min--snr_max
iterative_length = iterative_length+1
%------------(省略CP頭和訓練序列)上信道,乘性噪聲------------------------%

rec_code = zeros(Time_length,Rx_num,bits);
ch_mul_noice = zeros(Tx_num,Rx_num,bits);
ch_mul_noice = sqrt(0.5)*(randn(Tx_num,Rx_num,bits)+j*randn(Tx_num,Rx_num,bits));
for n = 1:Time_length 
for m = 1:bits
rec_code(n,:,m) = stbc_code(n,1:Tx_num,m)*ch_mul_noice(1:Tx_num,1:Rx_num,m);
end
end

%---------------------信號功率統計,加高斯白噪聲----------------------------%
% 注意：這里只有0.75的信號平均功率!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
rec_code = rec_code*sqrt(4/3);
sig_pante_pow = 0;
for n = 1:Rx_num
Agwn = sqrt(10^(-snr_db_step/10)/2)*(randn(Time_length,1,bits)+j*randn(Time_length,1,bits));
rec_code(:,n,:) = rec_code(:,n,:) + Agwn; 
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%---------------------接收信號并進行stbc譯碼(定點仿真開始)-----------------%

rec_code = rec_code/sqrt(4/3);
rec_code = rec_code*ch_scale; % 接收信號（疊加高斯白噪聲后）
[max_rec_sig,max_index] = max(max(max(max(real(rec_code),imag(rec_code)))));
max_rec_sig = fix(max_rec_sig/(fix(max_rec_sig)+1)*Max_value)/Max_value*(fix(max_rec_sig)+1);
rec_code = rec_code/max_rec_sig; % 信號歸一化,歸一化后信號都為純小數
rec_code = fix(rec_code*Max_value);

[max_ch_est,max_index] = max(max(max(max(real(ch_mul_noice),imag(ch_mul_noice)))));
max_ch_est = fix(max_ch_est/(fix(max_ch_est)+1)*Max_value)/Max_value*(fix(max_ch_est)+1);
ch_mul_noice = ch_mul_noice/max_ch_est; % 信道估計值歸一化
ch_mul_noice = fix(ch_mul_noice*Max_value);

rec_code_r = real(rec_code);
rec_code_i = imag(rec_code);
ch_mul_noice_r = real(ch_mul_noice);
ch_mul_noice_i = imag(ch_mul_noice);
decode_out = zeros((Tx_num-1),bits); % 輸出星座點的值
demap_int = zeros(Tx_num,bits); % 輸出qam映射的下標

for n = 1:bits
decode_tempt_r = zeros(1,(Tx_num-1)); % 信號長度=Tx_num （4）-1
decode_tempt_i = zeros(1,(Tx_num-1));
decode_tempt = zeros(1,(Tx_num-1));
m = reshape(ch_mul_noice(:,:,n),1,16); % 16 = Tx_num*Rx_num=4*4 ------ ch_mul_noice
tempt_eff = scale_eff;
for r = 1:Rx_num 
mult_tempt(1,1) = rec_code_r(1,r,n)*ch_mul_noice_r(1,r,n);
mult_tempt(1,2) = rec_code_i(1,r,n)*ch_mul_noice_i(1,r,n);
mult_tempt(1,3) = rec_code_r(2,r,n)*ch_mul_noice_r(2,r,n);
mult_tempt(1,4) = rec_code_i(2,r,n)*ch_mul_noice_i(2,r,n);
mult_tempt(1,5) = -rec_code_r(4,r,n)*ch_mul_noice_r(4,r,n);
mult_tempt(1,6) = -rec_code_i(4,r,n)*ch_mul_noice_i(4,r,n);
mult_tempt(1,7) = -rec_code_r(3,r,n)*ch_mul_noice_r(3,r,n);
mult_tempt(1,8) = -rec_code_i(3,r,n)*ch_mul_noice_i(3,r,n);
for nnn = 1:8
if mult_tempt(1,nnn) >= Max_add_value 
mult_tempt(1,nnn) = Max_add_value-1;
mult_flow_num = mult_flow_num+1;
elseif mult_tempt(1,nnn) < -Max_add_value
mult_tempt(1,nnn) = -Max_add_value;
mult_flow_num = mult_flow_num+1;
end
end

sum_tempt = fix((mult_tempt(1,1)+mult_tempt(1,2)+mult_tempt(1,3)+mult_tempt(1,4)+... 
mult_tempt(1,5)+mult_tempt(1,6)+mult_tempt(1,7)+mult_tempt(1,8))/2^tempt_eff);
decode_tempt_r(1,1) = decode_tempt_r(1,1)+sum_tempt; 

mult_tempt(1,1) = -rec_code_r(1,r,n)*ch_mul_noice_i(1,r,n);
mult_tempt(1,2) = rec_code_i(1,r,n)*ch_mul_noice_r(1,r,n);
mult_tempt(1,3) = rec_code_r(2,r,n)*ch_mul_noice_i(2,r,n);
mult_tempt(1,4) = -rec_code_i(2,r,n)*ch_mul_noice_r(2,r,n);
mult_tempt(1,5) = -rec_code_r(4,r,n)*ch_mul_noice_i(3,r,n);
mult_tempt(1,6) = rec_code_i(4,r,n)*ch_mul_noice_r(3,r,n);
mult_tempt(1,7) = -rec_code_r(3,r,n)*ch_mul_noice_i(4,r,n);
mult_tempt(1,8) = rec_code_i(3,r,n)*ch_mul_noice_r(4,r,n); 
for nnn = 1:8
if mult_tempt(1,nnn) >= Max_add_value 
mult_tempt(1,nnn) = Max_add_value-1;
mult_flow_num = mult_flow_num+1;
elseif mult_tempt(1,nnn) < -Max_add_value
mult_tempt(1,nnn) = -Max_add_value;
mult_flow_num = mult_flow_num+1;
end
end
sum_tempt = (mult_tempt(1,1)+mult_tempt(1,2) +...
mult_tempt(1,3)+mult_tempt(1,4)+...
mult_tempt(1,5)+mult_tempt(1,6)+...
mult_tempt(1,7)+mult_tempt(1,8))/2^tempt_eff;
decode_tempt_i(1,1) = decode_tempt_i(1,1)+sum_tempt;

mult_tempt(1,1) = rec_code_r(1,r,n)*ch_mul_noice_r(2,r,n);
mult_tempt(1,2) = rec_code_i(1,r,n)*ch_mul_noice_i(2,r,n);
mult_tempt(1,3) = -rec_code_r(2,r,n)*ch_mul_noice_r(1,r,n);
mult_tempt(1,4) = -rec_code_i(2,r,n)*ch_mul_noice_i(1,r,n);
mult_tempt(1,5) = rec_code_r(4,r,n)*ch_mul_noice_r(3,r,n);
mult_tempt(1,6) = rec_code_i(4,r,n)*ch_mul_noice_i(3,r,n);
mult_tempt(1,7) = -rec_code_r(3,r,n)*ch_mul_noice_r(4,r,n);
mult_tempt(1,8) = -rec_code_i(3,r,n)*ch_mul_noice_i(4,r,n);
for nnn = 1:8
if mult_tempt(1,nnn) >= Max_add_value 
mult_tempt(1,nnn) = Max_add_value-1;
mult_flow_num = mult_flow_num+1;
elseif mult_tempt(1,nnn) < -Max_add_value
mult_tempt(1,nnn) = -Max_add_value;
mult_flow_num = mult_flow_num+1;
end
end
sum_tempt = fix((mult_tempt(1,1)+mult_tempt(1,2) +...
mult_tempt(1,3)+mult_tempt(1,4)+...
mult_tempt(1,5)+mult_tempt(1,6)+...
mult_tempt(1,7)+mult_tempt(1,8))/2^tempt_eff); 
decode_tempt_r(1,2) = decode_tempt_r(1,2)+sum_tempt;

mult_tempt(1,1) = -rec_code_r(1,r,n)*ch_mul_noice_i(2,r,n);
mult_tempt(1,2) = rec_code_i(1,r,n)*ch_mul_noice_r(2,r,n);
mult_tempt(1,3) = -rec_code_r(2,r,n)*ch_mul_noice_i(1,r,n);
mult_tempt(1,4) = rec_code_i(2,r,n)*ch_mul_noice_r(1,r,n);
mult_tempt(1,5) = rec_code_r(4,r,n)*ch_mul_noice_i(4,r,n);
mult_tempt(1,6) = -rec_code_i(4,r,n)*ch_mul_noice_r(4,r,n);
mult_tempt(1,7) = -rec_code_r(3,r,n)*ch_mul_noice_i(3,r,n);
mult_tempt(1,8) = rec_code_i(3,r,n)*ch_mul_noice_r(3,r,n);
for nnn = 1:8
if mult_tempt(1,nnn) >= Max_add_value 
mult_tempt(1,nnn) = Max_add_value-1;
mult_flow_num = mult_flow_num+1;
elseif mult_tempt(1,nnn) < -Max_add_value
mult_tempt(1,nnn) = -Max_add_value;
mult_flow_num = mult_flow_num+1;
end
end
sum_tempt = fix((mult_tempt(1,1)+mult_tempt(1,2) +...
mult_tempt(1,3)+mult_tempt(1,4)+...
mult_tempt(1,5)+mult_tempt(1,6)+...
mult_tempt(1,7)+mult_tempt(1,8))/2^tempt_eff);
decode_tempt_i(1,2) = decode_tempt_i(1,2)+sum_tempt;

mult_tempt(1,1) = rec_code_r(1,r,n)*ch_mul_noice_r(3,r,n);
mult_tempt(1,2) = rec_code_i(1,r,n)*ch_mul_noice_i(3,r,n);
mult_tempt(1,3) = rec_code_r(2,r,n)*ch_mul_noice_r(3,r,n);
mult_tempt(1,4) = rec_code_i(2,r,n)*ch_mul_noice_i(3,r,n);
mult_tempt(1,5) = rec_code_r(1,r,n)*ch_mul_noice_r(4,r,n);
mult_tempt(1,6) = rec_code_i(1,r,n)*ch_mul_noice_i(4,r,n);
mult_tempt(1,7) = -rec_code_r(2,r,n)*ch_mul_noice_r(4,r,n);
mult_tempt(1,8) = -rec_code_i(2,r,n)*ch_mul_noice_i(4,r,n);
mult_tempt(1,9) = rec_code_r(3,r,n)*ch_mul_noice_r(1,r,n);
mult_tempt(1,10) = rec_code_i(3,r,n)*ch_mul_noice_i(1,r,n);
mult_tempt(1,11) = rec_code_r(3,r,n)*ch_mul_noice_r(2,r,n);
mult_tempt(1,12) = rec_code_i(3,r,n)*ch_mul_noice_i(2,r,n);
mult_tempt(1,13) = rec_code_r(4,r,n)*ch_mul_noice_r(1,r,n);
mult_tempt(1,14) = rec_code_i(4,r,n)*ch_mul_noice_i(1,r,n);
mult_tempt(1,15) = -rec_code_r(4,r,n)*ch_mul_noice_r(2,r,n);
mult_tempt(1,16) = -rec_code_i(4,r,n)*ch_mul_noice_i(2,r,n);
for nnn = 1:16
if mult_tempt(1,nnn) >= Max_add_value 
mult_tempt(1,nnn) = Max_add_value-1;
mult_flow_num = mult_flow_num+1;
elseif mult_tempt(1,nnn) < -Max_add_value
mult_tempt(1,nnn) = -Max_add_value;
mult_flow_num = mult_flow_num+1;
end
end
sum_tempt = fix((mult_tempt(1,1)+mult_tempt(1,2) +...
mult_tempt(1,3)+mult_tempt(1,4)+...
mult_tempt(1,5)+mult_tempt(1,6)+...
mult_tempt(1,7)+mult_tempt(1,8)+...
mult_tempt(1,9)+mult_tempt(1,10) +...
mult_tempt(1,11)+mult_tempt(1,12)+...