gardner1.m

通信中位同步的gardner算法的matlab仿真

clear
clc
close all

%     **********************  加载波和位定时同步   ***************************************************


fs = 88;    % 整个系统的采样频率
dt = 1/fs;
f_i = 22;
pi2 = 2*pi;
% c1 = 1/256;
% c2 = 1/(2048+128);
% c3 = 1/(2048+128);
c1 = 1/128;
c2 = 1/(2048+128);
c3 = 1/2048;     % 针对-0.02情况。
barkercode_11 = [1 1 1 -1 -1 -1 1 -1 -1 1 -1];
pn_15 = [1 1 1 1 0 1 0 1 1 0 0 1 0 0 0];
pn_sequence = pn_15*2-1;
coeff_5 = firrcos(28,5.5,0.25,88,'rolloff','sqrt');
% coeff_44 =[4 1 -4 -4 2 12 14 -2 -30 -48 -24 48 152 248 284 248 152 48 -24 -48 -30 -2 14 12 2 -4 -4 1 4];
Match_barker = kron(barkercode_11,ones(1,2));
c_1 = inline('1/6*u^3-1/6*u');
c_2 = inline('-1/2*u^3+1/2*u^2+u');
c_3 = inline('1/2*u^3-u^2-1/2*u+1');
c_4 = inline('-1/6*u^3+1/2*u^2-1/3*u');

fid = fopen('rand_data.m','r');
source1 =  fread(fid)';
fclose(fid);
source = [source1 zeros(1,100)];
sig_i = [1 pn_15 source(1:2:end)];
sig_q = [1 pn_15 source(2:2:end)];
match_regist = zeros(2,size(Match_barker,2));
% match_regist1 = zeros(2,size(Match_barker,2));
num = [0.01 1];
den = [1 1.01 1];
[Af,Bf,Cf,Df] = tf2ss(num,den);

for i = 1:length(sig_i)
    if i == 1
        DQpsk(1,i) = 1;
        DQpsk(2,i) = 1;
    else
        DQpsk(1,i) = abs(xor(sig_i(i),sig_q(i))-1)*xor(sig_i(i),pre_i) + xor(sig_i(i),sig_q(i))*xor(sig_q(i),pre_q);
        DQpsk(2,i) = abs(xor(sig_i(i),sig_q(i))-1)*xor(sig_q(i),pre_q) + xor(sig_i(i),sig_q(i))*xor(sig_i(i),pre_i);
    end
    pre_i = DQpsk(1,i);
    pre_q = DQpsk(2,i);
end
DQpsk = DQpsk*2-1;

Expand_sig = kron(DQpsk,barkercode_11);         % spread
Ep_insert_sig = zeros(2,length(Expand_sig)*8);  % insert 8 zero
Ep_insert_sig(:,1:8:end) = Expand_sig;

Filter_sig(1,:) = conv(Ep_insert_sig(1,:),coeff_5);  % low_filter (fc = 11,fs =88);
Filter_sig(2,:) = conv(Ep_insert_sig(2,:),coeff_5);

t = [0:dt:(size(Filter_sig,2)-1)*dt];
Md_sig_c = Filter_sig(1,:).*cos(pi2*22*t);
Md_sig_s = Filter_sig(2,:).*sin(pi2*22*t);   
clear Ep_insert_sig;
clear Filter_sig;
clear Ep_insert_sig;


sig_stor_c = resample(Md_sig_c,8798,8800);
sig_stor_s = resample(Md_sig_s,8798,8800);
dt1 = 1/87.98;
t = [0:dt1:(length(sig_stor_s)-1)*dt1];

noice = randn(1,length(sig_stor_c));
step = 2;
snr_min =  -4;
snr_max =  4;

err_rate = [];
store = [];
sig_ave_pow = (norm(sig_stor_c)^2+norm(sig_stor_s)^2)/length(sig_stor_c);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%  接受端  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
for snr = snr_min:step:snr_max
    snr
    Md_sig = sig_stor_c + sig_stor_s;  % I
    Md_sig = Md_sig/sqrt(sig_ave_pow) +sqrt(10^(-snr/10))*noice*sqrt(0.5);
   
    %--------------高频部分略过,通过下载波,可以得到 Md_sig_c 和 Md_sig_s --------------------
    nco_out = 0;
    nco_nature = 22.01;
    i = 0;
    dc_rig = zeros(1,29);
    ds_rig = zeros(1,29);
    match_out = zeros(2,22);
    phase_err = 0;
    fx = zeros(2,2);
    %-----------------Timing_parament--------------------------------------------
    it_flag = 0;
    count = 0;
    count1= 0;
    u = 0;
    timing_err = 0;
    match_count = 0;
    tnco_out = 1;
    interp_regist = zeros(2,4);  % 滤波器的存储器，将4个点存储。如果得到插值命令就用这4个值插值
    timing_regist_i = zeros(2,2);
    timing_regist_q = zeros(2,2);
    %     Ts = 1/22; 
    xt1 = zeros(1,2);
    xt2 = zeros(1,2);
    xt3 = zeros(1,2);
    contr_word = [0.25,0];
    %--------------------------------------------------------------------------  
    count2  = 0;
    pn_regist = zeros(2,15);
    pn_mul_result = 0;
    flag =0;
    pn_mul1 = 0;
    xf1 = zeros(1,2);
    xf2 = zeros(1,2);
    xf3 = zeros(1,2);
    decode_out = [];
    show_flag = 1;
    p_u = 0;
%     t = [0:dt:(length(sig_stor_s)-1)*dt];
   for i = 1:length(sig_stor_c)               
        %---------------- down convert frequency-----------------------------------------------
        
        dc = Md_sig(i)*cos(pi2*nco_nature*t(i));
        ds = Md_sig(i)*sin(pi2*nco_nature*t(i));
        dc_rig(2:end) = dc_rig(1:end-1);
        ds_rig(2:end) = ds_rig(1:end-1);
        dc_rig(1) = dc;
        ds_rig(1) = ds;
        
        d_convert_c = dc_rig*coeff_5';
        d_convert_s = ds_rig*coeff_5';
        
        %------------------ Timing synchronization------------------------------------------------
        
        count = count+1;                                % 调整计数器，反映从滤波器出来的数据个数，（包括插值的数值）
        interp_regist(:,2:end) = interp_regist(:,1:end-1);  % 滤波器的存储器，将4个点存储。如果得到插值命令就用这4个值插值
        interp_regist(1,1) = d_convert_c;
        interp_regist(2,1) = d_convert_s;
        if it_flag == 1      
            if xt2(1)==0
                u=0;
            else
                if abs(p_u-0.25)<10^-10 | abs(p_u+0.25)<10^-10
                    u = 0;
                else
                    u = p_u/contr_word(1)-1; 
                end
            end 
            interp_val_i = c_1(u)*interp_regist(1,4)+c_2(u)*interp_regist(1,3)... % 得到nco的指令，要求插值
                +c_3(u)*interp_regist(1,2)+c_4(u)*interp_regist(1,1);
            interp_val_q = c_1(u)*interp_regist(2,4)+c_2(u)*interp_regist(2,3)... % 得到nco的指令，要求插值
                +c_3(u)*interp_regist(2,2)+c_4(u)*interp_regist(2,1);
        end
        
        if mod(count,4)==1                                % 数据滤波，将数据率下降到原来的1/4
            count1 = mod(count1+1,2);                     % count1,寻找准点采样值.              
            if count1 == 0
                timing_regist_i(1,2) = timing_regist_i(1,1);
                timing_regist_i(1,1) = d_convert_c;     
                timing_regist_q(1,2) = timing_regist_q(1,1);
                timing_regist_q(1,1) = d_convert_s;
                
                data_filr_i = timing_regist_i(1,1);
                data_filr_q = timing_regist_q(1,1);
                %             timing_err = (timing_regist_i(1,1))/(abs(timing_regist_i(1,2))+abs(timing_regist_i(1,1)));
            elseif count1 ==1  
                timing_regist_i(2,2) = timing_regist_i(2,1);
                timing_regist_i(2,1) = d_convert_c;
                timing_regist_q(2,2) = timing_regist_q(2,1);
                timing_regist_q(2,1) = d_convert_s;
                
                data_filr_i = timing_regist_i(2,1);
                data_filr_q = timing_regist_q(2,1);
                
                timing_err = (timing_regist_i(2,2)*(timing_regist_i(1,2)-timing_regist_i(1,1))+...
                    timing_regist_q(2,2)*(timing_regist_q(1,2)-timing_regist_q(1,1)));  % ? 这里应该除二吗？？？
                %                     timing_err = (timing_regist_i(2,2)*(timing_regist_i(1,2)-timing_regist_i(1,1))+...
                %                         timing_regist_q(2,2)*(timing_regist_q(1,2)-timing_regist_q(1,1))); 
                %             timing_err = (timing_regist_i(2,1))/(abs(timing_regist_i(2,2))+abs(timing_regist_i(2,1))); 
                
            end     
        end
        
        xt1(2) = xt1(1);
        xt1(1) = c3*timing_err;  % 相位误差
        xt2(2) = xt2(1);
        xt3(2) = xt3(1);
        contr_word(2)=contr_word(1);
        
        xt2(1) = ((c1+c2)*xt1(1)-c1*xt1(2));  % 环路滤波器的输出
        contr_word(1) = contr_word(2)-xt2(1); % 得到控制字。W(m)
%         xt3(1) = xt2(1)+xt3(2)-0.25;                     % nco的输出

% tempt_sign = sign(xt2(1));
% if tempt_sign<0
    tnco_out = tnco_out-(contr_word(1));  
    p_u = tnco_out;
% else
%     p_u = tnco_out;
%     tnco_out = tnco_out-(contr_word(1));  
% end
        %             store = [store nco_out];
        if tnco_out<-10^-10
            it_flag =1;
            tnco_out = mod(tnco_out,1);
            continue;
        end
        
        %-----------------------------------------------------------------------
        
        if it_flag ==1
            it_flag =0;
            match_count = match_count + 1;
            %---------------Match filter -----------------------------------------------------------
            
            match_regist(:,1:end-1) = match_regist(:,2:end);
            match_regist(1,end) = data_filr_i;
            match_regist(2,end) = data_filr_q;
            
            match_out(:,2:end) = match_out(:,1:end-1);
            match_out(:,1) = [match_regist(1,:)*Match_barker(1,:)';match_regist(2,:)*Match_barker(1,:)']; 
            if mod(match_count,22) == 0
                [value index] = max(abs(match_out(1,:)));
                match_max_c = match_out(1,index);
                [value index] = max(abs(match_out(2,:)));
                match_max_s = match_out(2,index);
                
                tempt1 = match_max_c/sqrt(match_max_c^2+match_max_s^2);
                tempt2 = match_max_s/sqrt(match_max_c^2+match_max_s^2);
                
                match_max_c = tempt1;
                match_max_s = tempt2;
                
                if count2  ==0
                    diterm_sig_c = 1;
                    diterm_sig_s = 1;
                    S_afc = 0;
                    phase_err = S_afc;
                    %                         decode_out = [1 1];
                    pre_crot = sqrt(2)/2;
                    pre_srot = 0;
                    count2  = count2 +1;
                end
                if count2 >=1
                    count2  = count2 +1;
                    if count2  ==20
                        count2  = count2 ;
                    end
                    diterm_sig_c = match_max_c * pre_crot + match_max_s * pre_srot;
                    diterm_sig_s = match_max_s * pre_crot - match_max_c * pre_srot;
                    diterm_sig_c = diterm_sig_c/sqrt(2);
                    diterm_sig_s = diterm_sig_s/sqrt(2);
                    
                    pre_crot = (match_max_c+match_max_s)/2;
                    pre_srot = (match_max_s-match_max_c)/2;                   
                    
                    if    diterm_sig_c>0 & diterm_sig_s>0
                        decode(1) = -1;
                        decode(2) = -1;
                    elseif diterm_sig_c<0 & diterm_sig_s>0
                        decode(1) = -1;
                        decode(2) = 1;
                    elseif diterm_sig_c<0 & diterm_sig_s<0
                        decode(1) = 1;
                        decode(2) = 1;
                    elseif diterm_sig_c>0 & diterm_sig_s<0
                        decode(1) = 1;
                        decode(2) = -1;
                    end
                    if flag ==1  % flag 为帧同步标志位。
                        if decode(1)==-1
                            decode(1) = 0;
                        end
                        if decode(2) ==-1
                            decode(2) = 0;
                        end                              
                        decode_out = [decode_out decode];
                    else
                        pn_regist(:,1:end-1) = pn_regist(:,2:end);
                        pn_regist(:,end) = [decode(1);decode(2)];
                        pn_mul_result = pn_regist(1,:)*pn_sequence'+pn_regist(2,:)*pn_sequence';
                    end
                    %----------------  err inspect---------------------------
                    S_afc = sign(diterm_sig_c)*diterm_sig_s - sign(diterm_sig_s)*diterm_sig_c;
                    phase_err = S_afc/(match_max_c^2+match_max_s^2)/2;
                end         
                
            end
            
            if snr<=-8
                if pn_mul_result >=7
                    flag = 1;
                    if flag == 1&show_flag==1
                        pn_mul_result
                        show_flag =0;
                    end
                end
            elseif snr<-4
                if pn_mul_result >=10
                    flag = 1;
                    if flag == 1&show_flag==1
                        pn_mul_result
                        show_flag =0;
                    end
                end
            elseif snr<0
                if pn_mul_result >=15
                    flag = 1;
                    if flag == 1&show_flag==1
                        pn_mul_result
                        show_flag =0;
                    end
                end
            else
                if pn_mul_result >=20
                    flag = 1;
                    if flag == 1&show_flag==1
                        pn_mul_result
                        show_flag =0;
                    end
                end
            end
        else
            continue;
        end
        if i<=88
            nco_out = 0;
        else
            %             fx(:,2) = fx(:,1) +1.*Af*fx(:,1)/88 + 1.*Bf*phase_err/88;
            %             nco_out = Cf*fx(:,1);
            %             fx(:,1) = fx(:,2);
            xf1(2) = xf1(1);
            xf1(1) = phase_err-xf3(1);
            xf2(2) = xf2(1);
            xf3(2) = xf3(1);
            xf2(1) = (xf2(2)+(c1+c2)*xf1(1)-c1*xf1(2));  % 环路滤波器的输出
            xf3(1) = xf2(1)+xf3(2);
            nco_out = xf3(1);
            %             store = [store nco_out];
            
        end       
    end
    err_num = 0;
    err_bit = zeros(1,length(source1));
    if length(decode_out)==length(source1)-2
        decode_out = [decode_out 0 0];
    end
    for i = 1:length(source1)
        if (decode_out(i)~= source1(i))
            err_num = err_num+1;
            err_bit(i) = 1;
        end
    end
    err_rate = [err_rate err_num/length(source1)];
    err_num/length(source1)
end
% for i= 1:392
%     ee = decode_out(end);
%     decode_out(2:end) = decode_out(1:end-1);
%     decode_out(1) = ee;
%     a = length(find(decode_out~=source));
% end
% b=min(a)
snr = snr_min:step:snr_max;
semilogy(snr,err_rate,'r-*');
err_rate
grid on