try_ml.m

it is a simulation about ML synchronization algorithm in OFDM systems,you can slao see a function pi

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=128; % 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(1,N_symbol)*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.25;

for k=1:length(sn)
rx(k)=sn(k)*exp(j*2*pi*e*k/N); 
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
obs_rx=obs_rx+add_noise(obs_rx,SNR); % add noise
%################## 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);%equation(2-11)
gamma1(i,theta)=-angle(gamma(theta))/(2*pi);%frequency offset
estimate_d(i)=find(lamda(i,:)==max(lamda(i,:)));
end %end loop_theta

end %end if

cnt= cnt+1; % increase counter

end % end while


% S/P for lamda
lamda_str=[];
gamma1_str=[];
for i=1:N_symbol
estimate_d(i)=find(lamda(i,:)==max(lamda(i,:)));    
lamda_str=[ lamda_str lamda(i,:)];
gamma1_str=[ gamma1_str gamma1(i,:)];%set of frequency offset 
end
d_mvalue=mvalue(estimate_d(i))
delay_mvalue=mvalue(delay)
[d_sigma2,d_sigma]=disper(estimate_d(i))
d_bias=estimate_d-delay


% for display
clf;
grid on
subplot(2,1,1)
plot(lamda_str)
xlabel('时间');ylabel('最大似然函数的输出');
%plot(abs(xn)*100) 
subplot(2,1,2)
plot(gamma1_str)
xlabel('时间');ylabel('频率偏移输出');
error=gamma1_str(find(lamda_str==max(lamda_str)))-e