%本程序仿真載波恢復(fù)環(huán)
clc;
clear;
lengthdata = 2000
datai=randsrc(1,lengthdata,[-1 1]);
dataq=randsrc(1,lengthdata,[-1 1]);
stuff=zeros(1,15);
for i=1:length(datai)
    x((i-1)*16+1)=datai(i);
    x=[x stuff];
end;
data_i=x;
for i=1:length(dataq)
    y((i-1)*16+1)=dataq(i);
    y=[y stuff];
end;
data_q=y;

%平方根升余弦濾波器的系數(shù)
R=1.0;
n_T=[-4 4];
rate=16;
T=1;
b = rcosfir(R,n_T,rate,T,'sqrt');
wave_i=filter2(b,data_i,'same');
wave_q=filter2(b,data_q,'same');
I=filter2(b,wave_i,'same');
Q=filter2(b,wave_q,'same');

theta =pi/3; 
f=-2.4e3;                                                    %頻率偏移
fs=48e3;                                                    %采樣時鐘
m = 1:16:lengthdata*16-15;                   %最佳判決點所在的位置
signal = I(m) + j*Q(m) + 0.2*(rand(1,lengthdata)+j*rand(1,lengthdata));                         
k=1:lengthdata;
carrier = exp(j*2*pi*k*f/fs+theta);
complex = signal .* carrier ;
cplx = round(complex/max(abs(complex)).*(2^9-1));

Tsin=1;
fsin=1;
tsin=Tsin/512:Tsin/256:Tsin-Tsin/512;
x=sin(2*pi*fsin.*tsin);
sin_8bit=round(x.*(2^7-1));
Tcos=1;
fcos=1;
tcos=Tcos/512:Tcos/256:Tcos-Tcos/512;
y=cos(2*pi*fcos.*tcos);
cos_8bit=round(y.*(2^7-1));

% initilize PLL loop
phi_hat(1)=0;
e(1)=0;
phd_output(1)=0;
nco(1)=0;
address(1) = 1;


% define loop filter parameters
ki=1/256;      % proportional constant
kp=1/16;      % Integrator constant

% PLL implementation
for n=2:length(cplx)
    nco(n)=conj(cos_8bit(address(n-1)) + j*sin_8bit(address(n-1))); % Compute NCO
    xxx(n) = real(cplx(n)*nco(n));
    yyy(n) = imag(cplx(n)*nco(n));
    phd_output(n) = rem(round((sign(xxx(n))*yyy(n) - sign(yyy(n))*xxx(n))/2^0),2^16);
    e(n) = rem(round(ki * phd_output(n) + e(n-1)),2^16);
    phi_hat(n) = rem(round(phi_hat(n-1)+ e(n) + kp * phd_output(n)),2^16); 
    s(n) = mod(round(phi_hat(n)/2^8),256);
    address(n) = s(n) + 1-abs(sign(s(n)));
end;
% *************************************************


% plot waveforms
index_start = 1;
index_stop=1200;
figure
subplot(311),plot(index_start:index_stop, phd_output(index_start:index_stop)),ylabel('Ph. Det.'),grid on;
subplot(312),plot(index_start:index_stop, phi_hat(index_start:index_stop)*180/pi),ylabel('Est. Phs.'),grid on;
subplot(313),plot(index_start:index_stop,address(index_start:index_stop)),ylabel('address');grid on;
figure
subplot(211),plot(index_start:index_stop,real(nco(index_start:index_stop)),index_start:index_stop,real(carrier(index_start:index_stop)).*(2^6-1));
ylabel('RE-PLL')
subplot(212),plot(index_start:index_stop, imag(nco(index_start:index_stop)),index_start:index_stop,imag(carrier(index_start:index_stop)).*(2^6-1)),grid on;
ylabel('IM-PLL')

% plot(result);