clear;
Rb=5000000/31;%信息數(shù)據(jù)速率
lc=31;%擴(kuò)頻碼的周期。擴(kuò)頻增益13dB
Rc=Rb*lc;%碼元速率

L=1;%每個(gè)碼元里sin信號(hào)的周期個(gè)數(shù)
fc=Rc*L; %載波頻率

Q=8;%每個(gè)碼元中采樣個(gè)數(shù)
fs=Rc*Q;%采樣頻率
ts=1/fs;%每個(gè)取樣點(diǎn)的時(shí)間間隔

data_number=20;%信息數(shù)據(jù)的個(gè)數(shù)
sample_number=data_number*lc*Q;%采樣得到的總的數(shù)據(jù)個(gè)數(shù)
t=0:ts:(sample_number-1)*ts;

fo=fc;%經(jīng)過帶通采樣后的中心頻率
interference_band_init=fo-100000;%第一個(gè)干擾的起始位置
interference_band_final=fo+100000;%第一個(gè)干擾的終止位置
%%取均勻分布的隨機(jī)數(shù)據(jù)100個(gè)，即信息數(shù)據(jù)
temp=rand(1,data_number);
for i=1:data_number
    if (temp(i)<0.5)
        data1(i)=-1;
        data((i-1)*lc+1:i*lc)=-1;
    else
        data1(i)=1;
        data((i-1)*lc+1:i*lc)=1;
    end
end
%%生成周期為31位的m序列
register_current_value=[-1,-1,-1,-1,1];%5級(jí)移位寄存器
feed_position=[1,4];%反饋抽頭為1和4
L=length(register_current_value);
for i=1:data_number*lc  
    feed_value=0;
    for j=1:length(feed_position)
        feed_value=feed_value+register_current_value(feed_position(j));%計(jì)算反饋值
    end
    if feed_value==0
        register_new_value(L)=-1;
    else
        register_new_value(L)=1;
    end
    register_new_value(1:L-1)=register_current_value(2:L);
    register_current_value=register_new_value;
    pn(i)=register_current_value(1);
end
% figure
% stem(pn)
%%信息數(shù)據(jù)和PN碼相乘
sig=data.*1;
%%已調(diào)信號(hào)trans_sig
for i=1:data_number*lc
    resig((i-1)*Q+1:i*Q)=sig(i);
end   
a=1;%載波的幅度值
input_sig=a*resig.*cos(2*pi*(fc)*t);%調(diào)制

figure;
plot(input_sig);

%%載波同步的初始值
Wn=0.5*fc/fs;%截止頻率
ic=0.707;
k=ic*fc;
step=5;%濾波器階數(shù)
[B,A1]= butter(step,Wn,'low');%巴特沃思低通濾波器
A=A1(2:step+1);%a0為1，故后面不需要相乘
data_i_A=zeros(1,length(A));
data_i_B=zeros(1,length(B));
data_q_A=zeros(1,length(A));
data_q_B=zeros(1,length(B));
x_i=0;
x_q=0;
x_new=0;
y=0;
angle_value=0;%壓控振蕩器輸出的初始相位
%%載波同步的恢復(fù)
%科斯塔斯環(huán)
for i=1:length(input_sig)
    output_sig(i)=sin(2*pi*(fc-2000)*t(i)+angle_value);
    %分成兩路，一路乘以壓控振蕩器輸出的sin，另一路乘以壓控振蕩器輸出的cos
    e_i=input_sig(i)*cos(2*pi*(fc-2000)*t(i)+angle_value);%同相分量
    e_q=input_sig(i)*sin(2*pi*(fc-2000)*t(i)+angle_value);%正交分量
    
    nnn1(i)=e_i;
    nnn2(i)=e_q;
    %以下是經(jīng)過低通濾波器，濾除高頻分量
    data_i_B(1:length(B)-1)=data_i_B(2:length(B));%調(diào)整系數(shù)
    data_i_A(1:length(A)-1)=data_i_A(2:length(A));
    data_i_B(length(B))=e_i;
    data_i_A(length(A))=x_i;
    x_i=sum(data_i_B(length(B):-1:1).*B)-sum(data_i_A(length(A):-1:1).*A);%差分方程，實(shí)現(xiàn)低通濾波
    
    nn1(i)=x_i;
    
    data_q_B(1:length(B)-1)=data_q_B(2:length(B));%調(diào)整系數(shù)
    data_q_A(1:length(A)-1)=data_q_A(2:length(A));
    data_q_B(length(B))=e_q;
    data_q_A(length(A))=x_q;
    x_q=sum(data_q_B(length(B):-1:1).*B)-sum(data_q_A(length(A):-1:1).*A);%差分方程，實(shí)現(xiàn)低通濾波
       
    nn2(i)=x_q;
    %到此，完成濾波
    
    x_old=x_new;
    x_new=x_i*x_q;%兩路信號(hào)相乘，完成鑒相，得到相位差的正弦值x_new
    
    y=ts/2*(x_new+x_old)+y;%通過環(huán)路濾波器，得到控制信號(hào)y
    %壓控振蕩器
    angle_value=y*k;
    jieguo(i)=angle_value+pi/2;
end

figure
plot(jieguo)
figure
plot(input_sig)
figure
plot(output_sig)
figure
plot(1:length(t),input_sig,'r',1:length(t),output_sig,'b',1:length(t),jieguo,'g')