clear;
Rb=5000000/31;%信息數據速率
lc=31;%一個信息數據擴頻成31個碼元。擴頻增益13dB
Rc=Rb*lc;%碼元速率

L=13;%每個碼元里sin信號的周期個數
fc=Rc*L; %載波頻率

Q=4;%每個碼元中采樣個數
fs=Rc*Q;%采樣頻率
ts=1/fs;%每個取樣點的時間間隔

data_number=30;%信息數據的個數
sample_number=data_number*lc*Q;%采樣得到的總的數據個數
t=0:ts:(sample_number-1)*ts;

%%取均勻分布的隨機數據100個，即信息數據
temp=rand(1,data_number);
for i=1:data_number
    if (temp(i)<0.5)
        data1(i)=-1;
        data((i-1)*lc*Q+1:i*lc*Q)=-1;
    else
        data1(i)=1;
        data((i-1)*lc*Q+1:i*lc*Q)=1;
    end
end
a=1;%載波的幅度值
mod_sig=a*data.*cos(2*pi*fc*t);%調制，產生bpsk信號
%%位同步的初始值
[lowpass2_B,lowpass2_A1]=cheby2(3,10,0.1);%切比雪夫2型低通濾波器
lowpass2_A=lowpass2_A1(2:length(lowpass2_A1));%a0為1，所以乘的時候從a1開始乘即可
data_lowpass2_A=zeros(1,length(lowpass2_A));%經過低通濾波器后的數據
data_lowpass2_B=zeros(1,length(lowpass2_B));
output_lowpass2=0;

m=Q*lc/2;
N1=2;
N2=2; 
M1=3;

before_sum=0;
after_sum=0;
before_and_after_sum=0;
bit_before_division=0;
before_flag=0;
after_flag=0;

bit2_new=1;
bit2_old=0;
bit1_new=0;
bit1_old=1;

in_phase_integral=0;%同相積分值
middle_phase_integral=0;%中相積分值

in_phase_hold_new=0;
middle_phase_hold_new=0;

integral=0;
jjj=1;

pp=3;
e_q=mod_sig.*cos(2*pi*fc*t);%乘以相干載波
e_q=e_q(pp:length(e_q));%去掉前面2點
%畫圖
figure;
subplot(2,1,1);
plot(e_q);
subplot(2,1,2);
plot(filter(lowpass2_B,lowpass2_A1,e_q));
%%%%%%%%%%%%%%%%
for i=1:length(mod_sig)-11; %%位同步信號的恢復    
    
    %通過切比雪夫低通濾波器，濾除高頻分量
    data_lowpass2_B(1:length(lowpass2_B)-1)=data_lowpass2_B(2:length(lowpass2_B));%去掉前面1點
    data_lowpass2_A(1:length(lowpass2_A)-1)=data_lowpass2_A(2:length(lowpass2_A));%去掉前面1點
    data_lowpass2_B(length(lowpass2_B))=e_q(i);
    data_lowpass2_A(length(lowpass2_A))=output_lowpass2;
    output_lowpass2=sum(data_lowpass2_B(length(lowpass2_B):-1:1).*lowpass2_B)-sum(data_lowpass2_A(length(lowpass2_A):-1:1).*lowpass2_A);%差分方程的實現

    shuchu(i)=output_lowpass2;
    %判決，得到0、1方波序列
    if output_lowpass2>0
        output_lowpass2_reshape=1;
    else
        output_lowpass2_reshape=0;
    end   
    
    bit_input(i)=output_lowpass2_reshape;

    in_phase_hold_old=in_phase_hold_new;
    if bit2_new-bit2_old==1
        if in_phase_integral>0
            in_phase_hold_new=1;
        else
            in_phase_hold_new=0;
        end
        in_phase_integral=0;
    end
        
    if bit1_new-bit1_old==1
        if middle_phase_integral>0
            middle_phase_hold_new=1;
        else
            middle_phase_hold_new=0;
        end
        middle_phase_integral=0;
    end
    
    shuchu1(i)=bit1_new;
    shuchu2(i)=bit2_new;   
    
    bit1_old=bit1_new;
    bit2_old=bit2_new;
    
    in_phase_integral=in_phase_integral+bit_input(i);
    middle_phase_integral=middle_phase_integral+bit_input(i);
    
    in_phase_hold_diff=in_phase_hold_new-in_phase_hold_old;
    middle_phase_hold_xor_in=xor(in_phase_hold_new,middle_phase_hold_new);
    
    before=in_phase_hold_diff*middle_phase_hold_xor_in;
    after=in_phase_hold_diff*not(middle_phase_hold_xor_in);
    
    if before==1
        before_sum=before_sum+1;
        before_and_after_sum=before_and_after_sum+1;
    end
    if after==1
        after_sum=after_sum+1;
        before_and_after_sum=before_and_after_sum+1;
    end
    if before_sum==N1
        before_sum=0;
        after_sum=0;
        before_and_after_sum=0;
        bit_before_division=bit_before_division-1;%調整長度
    end
    if after_sum==N1
        before_sum=0;
        after_sum=0;
        before_and_after_sum=0;
        bit_before_division=bit_before_division+1;%調整長度
    end
    if before_and_after_sum==M1
        before_sum=0;
        after_sum=0;
        before_and_after_sum=0;
    end
    
    bit_before_division=bit_before_division+1;
    jlz(i)=bit_before_division;
    
    if bit_before_division>=m %分頻
        bit1_new=~bit1_old;
        bit2_new=~bit2_old;
        bit_before_division=rem(bit_before_division,m);
    else
        bit1_new=bit1_old;
        bit2_new=bit2_old;
    end
    bit_syn(i)=bit1_new-bit1_old;%找下跳邊沿
    if bit_syn(i)==-1
        if integral>Q*lc/2
            output_data(jjj)=1;
        else
            output_data(jjj)=-1;
        end
        integral=0;
        jjj=jjj+1;
    end
    integral=integral+output_lowpass2_reshape;       
end
%畫圖
figure
subplot(2,1,1)
stem(data(pp:length(data)))
subplot(2,1,2)
stem(shuchu2)