function [estNumerlogy,frmError] = STS_demod(rx_STF)
M=32;Q=16;OVR=1;
load STS_M32Q16.mat theta_mod v_mod STS_freq

STS_NUM = 9;
RX_re = reshape(rx_STF(17:end),[],STS_NUM);
STS_f = fftshift(fft(RX_re,16,1));
STS_freq_sk = STS_f([3:8,10:15],:);
theta_est = mod(angle(STS_freq_sk(2:end,:)./STS_freq_sk(1:end-1,:)),2*pi);

mmse_est = zeros(Q,STS_NUM);
for iii = 1:STS_NUM   
    for ii = 1:Q        
        mmse_tmp = abs(theta_est(:,iii)-theta_mod(ii,:).');   
        mmse_est(ii,iii) = norm(min([mmse_tmp,2*pi - mmse_tmp],[],2))^2;  
    end
end
mmse_est2 = sum(mmse_est,2);
v_est_no = find(mmse_est2==min(mmse_est2));
v_est = exp(j*v_mod(v_est_no(1)));
%%------demodulate phi 
STS_freq_sk = STS_freq_sk./abs(STS_freq_sk);        
% [~,STS_ref] = preamble_modu1(theta_mod(v_est_no(1),:), 0, OVR);
STS_ref = STS_freq(:,v_est_no(1),1);
STS_ref = STS_ref./abs(STS_ref);
phi_est_angle = mod(angle(sum(sum((STS_freq_sk(:,1:STS_NUM)./STS_ref)))),2*pi);
phi_est = exp(j*phi_est_angle);

% demap v and phi to Bits using pskmod 
v_data_est = pskdemod(v_est,Q,0,'gray');
dePreambleQ = de2bi(v_data_est,log2(Q));

phi_data_est = pskdemod(phi_est,M,0,'gray');
dePreambleM = de2bi(phi_data_est,log2(M));        

% CRC
poly = 'z4+z3+z2+z+1';
crcdetector = comm.CRCDetector(poly);
codeword = [dePreambleM,dePreambleQ];
[msg,frmError] = crcdetector(codeword.');
estNumerlogy = bi2de(msg.')+1;

end

%%
function [STS_time,STS_freq] = preamble_modu1(theta, phi, OVR)
% N_FFT =64;
s_k(1) = 1*exp(j*phi);
s_k(1) = sqrt(2)/abs(s_k(1))*s_k(1);
for i = 1:11
    s_k(i+1) = s_k(i)*exp(j*theta(i));
end
STS_LOC = [9,13,17,21,25,29,37,41,45,49,53,57];
Short_preamble_slot_Frequency = zeros(1,64); % [1x64]
Short_preamble_slot_Frequency(STS_LOC) = s_k;
STS_freq = ifftshift(Short_preamble_slot_Frequency);
STS = ifft([STS_freq(1:32),zeros(1,64*(OVR-1)),STS_freq(33:64)])*sqrt(64);
STS = sqrt(1/mean(abs(STS).^2))*STS;

STS_f = fft(STS);
STS_f = [STS_f(1:32),STS_f((end-31:end))];
STS_freq = fftshift(STS_f);

STS_time = STS(1:16*OVR);
STS_freq = STS_freq(STS_LOC);
% STS_dpsk = mod(angle(STS_freq((2:end))./STS_freq((1:end-1))),2*pi);
% STS_dpsk(theta==0) = min([STS_dpsk(theta==0);2*pi - STS_dpsk(theta==0)],[],1);
end