clear all; clc; close all;

SNR=-5:10;
N_FFT=1024;
NGI=1/4;
N_used=N_FFT/2;
N_GI=round(NGI*N_FFT);
Nofdm=N_FFT+N_GI;
Nsym = 500;

nSTOs=50;
CFO=0.4;

%% Формирование тренировочных символов
B=2;
reS=(-1).^round(rand(1,round(N_used/B)));
imS=1i.*(-1).^round(rand(1,round(N_used/B)));
s1=reS+imS;
s2=fliplr(s1);
tx_signal_noGI=[s1 s2 conj(s1) conj(s2)];
tx_signal = zeros(1,Nsym*N_FFT);
tx_signal_GI=[tx_signal_noGI(1,N_FFT-N_GI+1:N_FFT) tx_signal_noGI];
tx_signal=repmat(tx_signal_GI,1,Nsym);

%% Добавление смещения частоты и времени
if nSTOs>=0, y_STO=[tx_signal(1,nSTOs+1:end) zeros(1,nSTOs)];
else
y_STO=[zeros(1,-nSTOs) tx_signal(1,1:end+nSTOs)];
end
nn=0:length(y_STO)-1; y_CFO_STO = y_STO.*exp(1i*2*pi*CFO*nn/N_FFT);

%%
buff = [];
for k = 1:length(SNR)
recvd_signal=awgn(y_CFO_STO, SNR(k), 'measured', 'dB');
for n = 1:Nsym
%%
Len_all=N_FFT;
recvd_signal_zeropad=[zeros(1,Len_all) recvd_signal(1,(n-1)*Nofdm+1:n*Nofdm) zeros(1,Len_all) zeros(1,Len_all)];
for d=1:length(recvd_signal_zeropad)-Len_all-1
P_Park(d)=sum(recvd_signal_zeropad(d+(1:N_FFT/2)).*fliplr(recvd_signal_zeropad(d+N_FFT/2+(1:N_FFT/2)))); %
E_Park(d)=sum(abs(recvd_signal_zeropad(d+(1:N_FFT))).^2);
Q_Park(d)=sum(conj(recvd_signal_zeropad(d+(1:N_FFT/2))).*fliplr(recvd_signal_zeropad(d+N_FFT/2+(1:N_FFT/2))));
end
Park_metric=(abs(P_Park).^2)./(E_Park).^2;
buff = [buff Park_metric];
plot(Park_metric/max(Park_metric))

[~, maxipos1]=max(Park_metric-N_GI);
t_est_Park(1,n)=Nofdm - maxipos1;
FreqOffset(1,n) = angle(Q_Park(maxipos1-N_GI+1))/pi;
end
mse_CFO(k) = sqrt(sum(abs(FreqOffset - CFO).^2)/Nsym)
mse_STO(k) = sqrt(sum(abs(t_est_Park - nSTOs).^2)/Nsym)
end

clear all; clc; close all;

flTimeShift = 1; % 1 - с уходом времени; 0 - без ухода времени
flChann = 0; % 1 - с каналом; 0 - без канала

%choose SNR level in dB
SNR=-5:10;
%choose FFT length
N_FFT=1024;
%choose guard interval length as a percentage of N_FFT
NGI=1/4;
%choose number of used carriers
N_used=N_FFT/2;
N_GI=round(NGI*N_FFT);
Nofdm=N_FFT+N_GI;
Nsym = 500;
%frequency offset
nSTOs=50;
CFO=0.4;

Fs = 11.2e6; % Частота дискретизации сигнала
Tofdm = Nofdm/Fs; % Длительность OFDM сигнала в сек
ppm = 50*10^-6; % Стабильность тактового гененратора
Nppm = Fs*ppm; % Уход в отсчетах в сек

%% Канал 1
% PathDelays = [0 50e-9 110e-9 170e-9 290e-9 310e-9]; % задержки лучей
% U = [0 -3 -10 -18 -26 -32];
% Fd= [4.6 23 139]; % Доплеровское смещение в Гц
%% Канал 2
% PathDelays = [0 110e-9 190e-9 410e-9]; % задержки лучей
% U = [0 -9.7 -19.2 -22.8];
% Fd= [4.6 23 139]; % Доплеровское смещение в Гц
%% Канал 3
PathDelays = [0 310e-9 710e-9 1090e-9 1730e-9 2510e-9]; % задержки лучей
U = [0 -1.0 -9.0 -10.0 -15.0 -20.0];
Fd= [4.6 23 139]; % Доплеровское смещение в Гц

%% preamble Minn
B=2;
reS=(-1).^round(rand(1,round(N_used/B)));
imS=1i.*(-1).^round(rand(1,round(N_used/B)));
s1=reS+imS;
s2=fliplr(s1);
tx_signal_noGI=[s1 s2 conj(s1) conj(s2)];
tx_signal = zeros(1,Nsym*N_FFT);
tx_signal_GI=[tx_signal_noGI(1,N_FFT-N_GI+1:N_FFT) tx_signal_noGI];
tx_signal=repmat(tx_signal_GI,1,Nsym);

%% Модель канала
OutChan = complex(zeros(Nsym*N_FFT,1),zeros(Nsym*N_FFT,1));
rand('seed',0);
if flChann == 1
tic
%T = [0 PathDelays(OptChan)];
%U = [0 0];
T = PathDelays;
chan = rayleighchan(1/Fs,Fd(3)/2,T,U);
chan.DopplerSpectrum=doppler.ajakes; %flat ajakes
chan.ResetBeforeFiltering=0;
OutChan=filter(chan,tx_signal);
t = toc;
fprintf('Сигнал искажен каналом...\n Время выполнения %4.2f секунд\n', t);
else
OutChan = tx_signal;
end

%% Добавление смещение частоты и времени
tic;
delay(1,1) = Nppm*Tofdm;
L = Nofdm;
for n = 2:Nsym
delay(1,n) = delay(1,n-1) + Nppm*Tofdm;
int_delay = fix(delay);
fract_delay = delay - int_delay;
end

MaxDelay = fix(max(int_delay));

OutChan(1,end:end+MaxDelay) = 0;
buff_delay = complex(zeros(1, 2*L), zeros(1, 2*L));
if flTimeShift == 1
for n = 1:Nsym
buff_delay = circshift(buff_delay, [0 -(L+int_delay(1,n))]);
buff_delay(1,L-int_delay(1,n)+1:end) = OutChan(1,(n-1)*L+1:n*L+int_delay(1,n));

re_buff_delay = sig_delay(real(buff_delay),1/Fs,fract_delay(1,n)/Fs);
im_buff_delay = sig_delay(imag(buff_delay),1/Fs,fract_delay(1,n)/Fs);

buff_delay = re_buff_delay + 1i*im_buff_delay;

y_STO(1,(n-1)*L+1:n*L) = buff_delay(1,L+1:2*L);
end
else
y_STO = OutChan;
end

if nSTOs>=0, y_STO=[y_STO(1,nSTOs+1:end) zeros(1,nSTOs)];
else
y_STO=[zeros(1,-nSTOs) y_STO(1,1:end+nSTOs)];
end

nn=0:length(y_STO)-1; y_CFO_STO = y_STO.*exp(1i*2*pi*CFO*nn/N_FFT);
t = toc;
fprintf('Осуществлен сдвиг времени и частоты сигнала...\n Время выполнения %4.2f секунд\n', t);

%%
buff = [];
for k = 1:length(SNR)
recvd_signal=awgn(y_CFO_STO, SNR(k), 'measured', 'dB');

for n = 1:Nsym
%%
Len_all=N_FFT; %
recvd_signal_zeropad=[zeros(1,Len_all) recvd_signal(1,(n-1)*Nofdm+1:n*Nofdm) zeros(1,Len_all) zeros(1,Len_all)];
for d=1:length(recvd_signal_zeropad)-Len_all-1
P_Park(d)=sum(recvd_signal_zeropad(d+(1:N_FFT/2)).*fliplr(recvd_signal_zeropad(d+N_FFT/2+(1:N_FFT/2)))); %fliplr
E_Park(d)=sum(abs(recvd_signal_zeropad(d+(1:N_FFT))).^2);
Q_Park(d)=sum(conj(recvd_signal_zeropad(d+(1:N_FFT/2))).*fliplr(recvd_signal_zeropad(d+N_FFT/2+(1:N_FFT/2))));
end
Park_metric=(abs(P_Park).^2)./(E_Park).^2;
buff = [buff Park_metric];
plot(Park_metric/max(Park_metric))

[~, maxipos1]=max(Park_metric-N_GI);
t_est_Park(1,n)=Nofdm - maxipos1;
FreqOffset(1,n) = angle(Q_Park(maxipos1-N_GI))/pi;
end

mse_CFO(k) = sqrt(sum(abs(FreqOffset - CFO).^2)/Nsym)
if flTimeShift == 1
mse_STO(k) = sqrt(sum(abs(diff(t_est_Park - nSTOs) - diff(int_delay(1,)).^2)/Nsym) %
else
mse_STO(k) = sqrt(sum(abs(t_est_Park - nSTOs).^2)/Nsym)
end

end

function sig_d = sig_delay(sig,T,dt)
% Задержка сигнала по времени
% sig - отсчеты сигнала
%T - шаг дискретизации
%dt - задержка по времени
%sig_d - задержанный сигнал
if dt==0
sig_d=sig;
return
end

N=numel(sig);
sig_d(1:N)=0;
a=fix(dt/T);
b=dt/T-a;

x=-b;
y0=0;y1=0;y2=0;y3=0;

for m=1:N
if m-a-2>0
y0=sig(m-a-2);
end
if m-a-1>0
y1=sig(m-a-1);
end
if m-a>0
y2=sig(m-a);
end
if (m-a+1>0)&&(m-a+1<=N)
y3=sig(m-a+1);
end
a3=(y3-y0)/6+(y1-y2)/2;
a1=(y3-y1)/2-a3;
a2=y3-y2-a1-a3;
a0=y2;

sig_d(m)=((a3*x+a2)*x+a1)*x+a0;
end
end