none_timing_resample.m

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

clear
clc
close all

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


fs = 88;    % 整个系统的采样频率
dt = 1/fs;
f_i = 22;
pi2 = 2*pi;
c1 = 1;
c2 = 1/32;
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_11 = firrcos(28,11,0.25,88,'rolloff','sqrt');
coeff_22 = firrcos(28,22,0.25,88,'rolloff','sqrt');
% coeff_22 =[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.1 1];
den = [1 1.5 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 = 44);
Filter_sig(2,:) = conv(Ep_insert_sig(2,:),coeff_5);

% ppm = 176*10^-2;
% dt1 = 1/(88+ppm);
t = [0:dt:(size(Filter_sig,2)-1)*dt];

Md_sig_c =  Filter_sig(1,:).*cos(pi2*22*t) + Filter_sig(2,:).*sin(pi2*22*t);
Md_sig_s = -Filter_sig(1,:).*sin(pi2*22*t) + Filter_sig(2,:).*cos(pi2*22*t);   % 与希尔伯特变换等效,
clear Ep_insert_sig;
clear Filter_sig;
clear Ep_insert_sig;

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

noice = randn(1,length(sig_stor_c));
% noice = randn(1,length(Md_sig_c));
step = 2;
snr_min = -4;
snr_max = 4;
sig_ave_pow = (norm(Md_sig_c)^2+norm(Md_sig_s)^2)/(length(Md_sig_s));
% sig_stor_c = Md_sig_c;
% sig_stor_s = Md_sig_s;
err_rate = [];
store = [];


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%  接受端  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
for snr = snr_min:step:snr_max
    snr
    Md_sig_c = sig_stor_c./sqrt(sig_ave_pow) + sqrt(10^(-snr/10))*sqrt(0.5)*noice;
    Md_sig_s = sig_stor_s./sqrt(sig_ave_pow) + sqrt(10^(-snr/10))*sqrt(0.5)*noice;
% Md_sig_c = sig_stor_c/sqrt(sig_ave_pow) ;
% Md_sig_s = sig_stor_s/sqrt(sig_ave_pow);
    Md_sig = Md_sig_c + Md_sig_s;
    %--------------高频部分略过,通过下载波,可以得到 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);
    num = length(Md_sig_c);
    match_out = zeros(2,22);
    phase_err = 0;
    fx = zeros(2,2);
    
    count1 = 0;
    u = 0;
    timing_count = 0;
    tnco_out = 1;
    %     Ts = 1/22;
    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 = 0;
    while num>0
        num = num-1;       
        i = i+1;        
        %---------------- down convert frequency-----------------------------------------------
        
        dc = Md_sig_c(i)*cos(pi2*nco_nature*t(i)) - Md_sig_s(i)*sin(pi2*nco_nature*t(i));
        ds = Md_sig_c(i)*sin(pi2*nco_nature*t(i)) + Md_sig_s(i)*cos(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';
        
        if mod(i,4)==1
            timing_count = timing_count + 1;
            %------------------ Timing synchronization------------------------------------------------
            
            %---------------Match filter -----------------------------------------------------------
            
            match_regist(:,1:end-1) = match_regist(:,2:end);
            match_regist(1,end) = d_convert_c;
            match_regist(2,end) = d_convert_s;
            
            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(timing_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 count1 ==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;
                    count1 = count1+1;
                end
                if count1>=1
                    count1 = count1+1;
                    if count1 ==20
                        count1 = count1;
                    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
                        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<=-10
                if pn_mul_result >=10
                    flag = 1;
                    if flag == 1&show_flag==1
                        pn_mul_result
                        show_flag =0;
                    end
                end
            elseif snr<-6
                if pn_mul_result >=15
                    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 >=10
                    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)];
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,'b-o');
err_rate
grid on