clear all;
clc;

N_symbol=10;

N = 1024; %the number of carrier
%N = 256; %the number of carrier

%T=10^(-6); % Sampling time : 1us

L=512; % the number of prefix
SNR=15; % SNR of signal 

rx_pre2 = zeros(1,N+L); %buffer for symber i-1
rx_pre1 = zeros(1,N+L); %buffer for symbol i
rx = zeros(1,N+L); %buffer for symbol i+1

ro = SNR/(SNR+1);
%e=0.25;
%abs(e)=0;

delay=floor(rand*N); %generate delay

cnt = 0;
while cnt <= N_symbol+1

% Generate Signal
input_stream = (sign(randn(1,N*2))+1)/2;

% modulation to QPSK
qpsk_str = QPSK(input_stream); % QPSK stream : 1024

% inverse fft : from frequency domain to time domain
xn=ifft(qpsk_str,N);

% add_prefix
len_xn = length(xn); % for checking : len_xn=N
sn=[xn(len_xn-L+1:len_xn) xn];

rx_pre2=rx_pre1;
rx_pre1=rx;
%############################ AWGN Channel #############################
e=0;
for theta=1:N
for k=theta+1:length(sn)
rx(k)=sn(k-theta)*exp(sqrt(-1)*2*pi*e*k/N);
end
end
rx= rx+add_noise(sn,SNR); % add noise
%rx(delay:length(rx))=rx;
%rx(1:delay-1)=0;
obs_rx=[rx_pre2(N+L-delay+1:N+L) rx_pre1 rx(1:N-delay)]; % observe 2N+L

%################## Calculate gamma(m) & PI(m) ###################

if cnt > 1

i=cnt-1;
gamma =zeros(1,N);
pii =zeros(1,N);

for theta=1:N

%Equation 2-(6),2-(7)
for k=theta:theta+L-1
gamma(theta)= gamma(theta)+obs_rx(k)*conj(obs_rx(k+N));
pii(theta)=pii(theta)+0.5*( abs(obs_rx(k))^2 + abs(obs_rx(k+N))^2 );
end

%equation 2-(5)
%lamda(i,theta)=abs(gamma(theta))*cos( 2*pi*e+angle(gamma(theta)) ) - ro*pii(theta);
lamda(i,theta)=abs(gamma(theta)) - ro*pii(theta);

end %end loop_theta

end %end if

cnt= cnt+1; % increase counter

end % end while

% S/P for lamda
lamda_str=[];
for i=1:N_symbol
lamda_str=[ lamda_str lamda(i,：];
end

%off=find(lamda_str==max(lamda_str))-1;
%xnn=[];
%for i=1:N_symbol
% xnn=[xnn xn(1024*(i-1)+1:1024)];
%end
%err=xnn(1:length(sn)-off)-rx(off+1:length(sn));
%te=sum(abs(err));
%ber=te/length(sn);

% for display
clf;
%axis([0 2000 -0.5 0.1]);
figure(1);
grid on
plot(lamda_str)
hold on 

function [rx]= add_noise(sn,SNR)
fb=10^5;
fs=10^6;
ebn0=10^(SNR/10);
eb=sum(sn.^2)./(length(sn)*fb);
n0=eb/ebn0;
pn=n0*fs/2;
n=sqrt(pn)*randn(1,length(sn));
rx=sn+n;

function [S] = QPSK(bits)
nb = length(bits); % number of bits
ns = ceil(nb/2); % number of symbols
b0 = zeros(1,ns*2); % zero padding
b0(1:nb) = bits;
j = sqrt(-1);
for s = 1 : ns

ba = b0(((s-1)*2)+1);
bb = b0(((s-1)*2)+2);
k = bb + ba*2;
p = ((pi/4)*(2*k-1))-pi;

Sc(s) = cos(p);
Ss(s) = sin(p);
S(s) = Sc(s) + j*Ss(s);

end