clear all; clc; close all;

SNR=50;
N_FFT=1024;
NGI=1/4;

N_used=N_FFT/2;
N_GI=round(NGI*N_FFT);
N=N_FFT+N_GI;

nSTOs=0;
CFO=0;

%% preamble Minn
B=4;
reS=(-1).^round(rand(1,round(N_used/B)));
imS=1i.*(-1).^round(rand(1,round(N_used/B)));
qam_data_tx= reS+imS;
s=ifft(qam_data_tx, N_FFT/B);
tx_signal_noGI = [s s -s -s];
%tx_signal=[tx_signal_noGI(:,N_FFT-N_GI+1:N_FFT) tx_signal_noGI];
tx_signal=tx_signal_noGI;
N=N_FFT+N_GI;
%% Time and frequency shift signal
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(j*2*pi*CFO*nn/N_FFT);

%% AWGN channel
recvd_signal=awgn(y_CFO_STO, SNR, 'measured', 'dB');
%% Minn
Len_all=N_FFT;
recvd_signal=[zeros(1,Len_all) recvd_signal zeros(1,Len_all)];

E_Minn=zeros(1,2*Len_all+1);
P_Minn=zeros(1,2*Len_all+1);
for d = 1:2*Len_all+1
for m = 0:1
rx_sgn_b = recvd_signal(d+m*Len_all/2:d+m*Len_all/2+Len_all/B-1);
rx_sgn_bp1 = recvd_signal(d+m*Len_all/2+Len_all/B:d+m*Len_all/2+Len_all/B+Len_all/B-1);
P_Minn(d)=P_Minn(d)+(sum(conj(rx_sgn_b).*rx_sgn_bp1));
E_Minn(d)=E_Minn(d)+sum(abs(rx_sgn_bp1).^2);
end
end
Minn_metric=(abs(P_Minn)).^2./(E_Minn).^2;

plot(Minn_metric);