% clc
clear all
close all

T=1;
L=3;
h=2/7;
m=8;            % 八進制
p = 2/h;        % 累計相位的狀態個數
OVS=40;
OVR=40;
EbN0 = 10;
SNR = EbN0 - 10 * log10(OVR) + 10 * log10(log2(m));
SNR = 10;
% SNR = 1;

block_num = 1;      %發送信息的楨數
block_size = 3;     %發送信息的楨大小   block_num*block_size為發送的信息bit數
frm_num = block_num;    
frm_size = block_size + 3;  %每楨加三個歸零碼

t=[0:T/OVS:L*T];
gt=( 1-cos( 2*pi*t/(L*T) ) )/(2*L*T);

for ii=1:L*OVS
    QT(ii)=sum(gt(1:ii).*(T/OVS));
end
rand('state',2);
%source_bit = 2*(randint(1,block_num*block_size,m)) - (m-1);     %生成發送信息
source_bit = [7,5,3];

% % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % %
% 對信息序列成楨，相當于每一條發的序列，在每一楨加三個symbol是狀態歸零
% % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % 
framed_bit = zeros(1,block_num*(block_size+3));
tail = [0,-(m-1),-(m-1)];
v = [];
for ii = 1:block_num
    v = source_bit( (ii-1)*block_size+1:ii*block_size );
    tail(1) = mod(p - mod(sum((v+m-1)/2),p),m);
    tail(1) = tail(1)*2 - (m-1);
    framed_bit((ii-1)*(block_size+3)+1:ii*(block_size+3)) = [source_bit((ii-1)*block_size+1:ii*block_size),tail];
end;


% % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % 
%                                       MODULATE
% % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % 
init_state = [0,0,0];
phas = zeros(1,length(framed_bit));
acc_ph = zeros(1,length(framed_bit)*OVS);
mod_bit = [-(m-1),-(m-1),framed_bit];	% init state = (0,-7,-7) 使第一楨從初始狀態（phase=0,ak2=-7,ak1=-7）開始調制

for jj = 3:length(mod_bit)
		kk = jj - 2;
	if(kk>1)
		phas(kk) = mod(phas(kk-1) + pi*h*mod_bit(jj-3),2*pi);
	end;
    acc_ph((kk-1)*OVS+1:kk*OVS) = mod(  2*pi*h*( mod_bit(jj-2)*QT(2*OVS+1:3*OVS) + mod_bit(jj-1)*QT(OVS+1:2*OVS) + mod_bit(jj)*QT(1:OVS) ) + phas(kk),2*pi );
end;
mod_sig = exp(j*acc_ph);
% % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % 
%     PLOT ACC_PH
% % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % 
% x = 1./OVS:1./OVS:length(mod_bit)-2;
% plot(x,acc_ph);
% ylim([0,2*pi]);

% % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % 
%											TRANSMIT
% % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % 

%為每一跳加上一個隨機相位
trans_sig1 = zeros(1,length(mod_sig));
trans_sig2 = zeros(1,length(mod_sig));
rand('state',3);
for ii = 1:frm_num;
	index1 = (ii-1)*frm_size*OVS + 1;
	index2 = ii*frm_size*OVS;
	trans_sig1(index1:index2) = mod_sig(index1:index2)*exp(j*1.9);
end;
trans_sig2 = mod_sig;    %   不加隨機相位


% % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % 
%											CHANNEL
% % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % %

noisy_sig1 = awgn(trans_sig1, SNR, 0, 1);
noisy_sig2 = awgn(trans_sig2, SNR, 0, 1);
noisy_sig1 = trans_sig1;    %沒有噪聲的情況  
noisy_sig2 = trans_sig2;
% Ir = It;
% Qr = Qt;


% % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % 
% %                                        DEMODULATE
% % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % 
rec_sig1 = noisy_sig1;
rec_sig2 = noisy_sig2;

v1 = rec_sig1(1:40)*mod_sig(1:40)';
abs(v1)
v2 = rec_sig1(1:40)*mod_sig(41:80)';
abs(v2)
v3 = rec_sig1(1:40)*mod_sig(81:120)';
abs(v3)

trans_ph1 = mod(real(log(trans_sig1)/j),2*pi);
trans_ph2 = mod(real(log(trans_sig2)/j),2*pi);
noisy_ph1 = mod(real(log(noisy_sig1)/j),2*pi);
noisy_ph2 = mod(real(log(noisy_sig2)/j),2*pi);
rec_ph1  = mod(real(log(rec_sig1)/j),2*pi);
rec_ph2  = mod(real(log(rec_sig2)/j),2*pi);

x = 1./OVS:1./OVS:length(mod_bit)-2;
figure(1);
subplot(4,1,1);plot(x,acc_ph);ylim([0,2*pi]);
subplot(4,1,2);plot(x,trans_ph1);ylim([0,2*pi]);
subplot(4,1,3);plot(x,noisy_ph1);ylim([0,2*pi]);
subplot(4,1,4);plot(x,rec_ph1);ylim([0,2*pi]);
% figure(2);
% subplot(4,1,1);plot(x,acc_ph);ylim([0,2*pi]);
% subplot(4,1,2);plot(x,trans_ph2);ylim([0,2*pi]);
% subplot(4,1,3);plot(x,noisy_ph2);ylim([0,2*pi]);
% subplot(4,1,4);plot(x,dmod_ph2);ylim([0,2*pi]);
hold on