freq_time_recovery.m

PHS多天线基带调制与解调MATLAB程序

function [align, freq] = Freq_Time_Recovery(signal, overSample, symRate);

begSym = 1;
numSym = length(signal) / overSample - 1;

apstart = (begSym-1)*overSample+1;
apend   = numSym*overSample + apstart -1;

z = signal(apstart:apend).';

k     =  length(z);

n= floor(k/overSample - 1);
next  =  (overSample+1):overSample*n+overSample;
curr  =  1:overSample*n; 


diffPhasor   =  z(next) .* conj(z(curr));


quad = 2*(imag(diffPhasor) < 0) + (real(diffPhasor) < 0) + 1;

rot = [5931641-5931641i -5931641-5931641i 5931641+5931641i -5931641+5931641i ];
%rot = [1-i -1-i 1+i -1+i ];



derotateDiffPhasor = floor(2 * (z(curr) .* conj(rot(quad))) / 2^24);
derotateDiffPhasor = reshape(derotateDiffPhasor, overSample, n);
zz = reshape(z(next),overSample, n);

zabsq = abs(z).^2;
 
zabsq = zabsq(1:overSample*floor(length(zabsq)/overSample));

idprod = kron(ones(length(zabsq)/overSample,1),eye(overSample));
zsum = (zabsq * idprod);
tone = exp(2*pi*j*[0:overSample-1]/overSample);

Tau_k = 1+angle(tone*zsum')/2/pi;
if(Tau_k > 1)
  Tau_k = Tau_k-1;
end

align = rem(round(Tau_k * overSample),overSample);
maxidx = align + 1;

align = align + 1;

derotateDiffPhasor = diffPhasor .* rot(quad);
derotateDiffPhasor = reshape(derotateDiffPhasor, overSample, n);
tm = derotateDiffPhasor(maxidx, :);

%ap = 2^23 - 1;
%an = 2^23 - 1;
ap = 1;
an = 1;

cntp = 0;
cntn = 0;
for i=1:length(tm)
   ti = tm(i);
   
if imag(ti) > 0
   
   old_ap = real(ap) * imag(ap);
   
	ap = ap * ti;
	apabs = max(abs([real(ap) imag(ap)]));
   if (apabs==0) clb = 23;
   else clb = floor(23 - log(apabs)/log(2));
   end
	ap = ap * 2^clb;
	new_ap = real(ap) * imag(ap);


	if old_ap*new_ap < 0, cntp = cntp + 1; end
    else


	old_an = real(an) * imag(an);


	an = an * ti;
	anabs = max(abs([real(an) imag(an)]));
   if (anabs ==0) clb = 23;
   else clb = floor(23 - log(anabs)/log(2));
   end
	an = an * 2^clb;
	new_an = real(an) * imag(an);

	if old_an*new_an < 0, cntn = cntn + 1; end
    end
end

rot2 = floor(2^23 * exp([0:2:6]*(-j*pi/4)) + 0.5);
ap = ap * rot2(rem(cntp, 4)+1);
an = an * conj(rot2(rem(cntn, 4)+1));
ap = ap/2;

apn = ap * an;
%t_apn = apn;
if imag(apn) > 0
    apn = apn * rot(1);
    cntp = cntp + 0.5;
else
    apn = apn * conj(rot(1));
    cntn = cntn + 0.5;
end

%k = 2;
%temp = 0;
%for i = 1:4
%    x = -j*sign(imag(apn))*pi/(2^k);
%    temp = temp + sign(imag(apn))*pi/(2^k);
%    x = exp(x);
%    apn = apn*x;
%    k = k+1;
%end
 
%t=1;
%for i=1:length(tm);
%    t = t*tm(i)/2^23;
%end

%temp = (temp + (cntp - cntn)*pi/2)/length(tm);
%temp = (temp + imag(apn)/real(apn) + (cntp - cntn)*pi/2)/length(tm);
temp = (angle(apn) + (cntp - cntn)*pi/2)/length(tm);
freq = temp*symRate/2/pi;
align = align - 1;