test_finally_edition.asv

此为频偏校正程序,通过TD_GG算法得到的频偏值还存在几赫兹的残余频偏

%%%%%%%%%%%%%%%%%%

clear;
global input
a4=[0 0 0;0 0 1;0 1 1;0 1 0;1 1 0;1 0 0;1 0 1;1 1 1];
%%%%%%%%%%%%%%%%%%%%%%%%%
A=[3*pi/4 pi pi/4 pi/2 -pi/4 -pi/2 -3*pi/4 0];%%%   A应该是正确的格雷码映射
%%%%%%%%%%%%%%%%%%%%%%%
M=8; 
T=1/270833; %%%%%%%%%%%%%%%%%%%%
Ts=T/4;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%
N=200;


snr_in_dB=10; 
Es=1;
snr=10^(snr_in_dB/10);
sgma=sqrt(Es/(sqrt(M))/snr/2);
%input=[1	1	1	0	1	1	1	0	1	0	1	1	1	0	0	1	1	1	1	0	1	0	1	1	1	0	0	0	0	0	1	1	0	1	0	1	1	0	1	0	1	1	1	1	0	1	1	0	0	0	1	0	1	0	0	1	0	0	0	1	1	0	1	1	1	1	1	1	1	0	0	1	1	0	1	1	0	0	0	0	0	0	1	0	1	0	1	0	1	0	0	1	0	0	1	1	0	1	0	0	0	0	1	1	0	1	0	0	0	1	1	0	0	1	1	1	0	0	1	1	0	0	0	1	1	0	1	0	1	0	0	1	0	1	1	1	1	1	0	0	0	1	0	0	0	1	1	0	0	0	0	0	1	0	0	0	1	1	0	1	1	1	1	1	0	1	0	0	0	0	1	1	1	1	1	1	1	0	1	0	1	0	1	1	1	1	1	1	1	0	1	0	0	1	0	1	1	1	1	0	1	1	0	1	1	0	1	0	1	0	1	0	0	0	0	0	0	0	0	0	1	0	0	0	0	0	1	0	1	0	1	0	1	0	1	1	1	0	1	0	1	1	0	1	0	1	1	1	1	0	0	0	0	1	0	0	0	0	1	1	1	0	1	0	0	0	1	1	0	1	1	1	1	0	1	0	0	0	1	0	0	1	1	0	1	1	1	1	1	1	1	1	0	1	0	0	1	1	0	0	1	1	0	1	0	0	0	1	1	1	1	0	0	1	0	1	1	0	0	1	0	1	1	1	1	1	0	1	1	0	1	0	0	1	0	1	1	1	1	1	1	0	0	1	1	0	1	1	1	1	1	0	0	1	0	1	0	0	0	1	0	1	1	1	0	0	1	0	0	0	1	1	1	1	0	1	1	0	1	1	0	1	1	1	0	0	1	1	0	0	0	0	1	1	1	1	0	1	1	1	0	0	0	1	1	0	1	0	0	0	1	1	1	1	1	1	0	0	0	0	0	1	0	0	0	1	1	0	0	1	1	1	0	1	1	0	1	0	0	0	1	1	0	0	1	0	1	1	0	1	1	0	1	0	0	0	0	1	1	1	0	1	1	0	1	1	1	1	0	0	0	1	1	0	1	1	1	1	1	0	0	0	0	1	0	0	1	1	0	0	1	0	1	0	0	0	0	1	0	1	1	1	1	1	1	0	0	1	1	1	0	0	1	1	0	1	0	0	1	0	1	0	0	0	1	1	1	1	1	1	0	1	0	0	0	0	1	0	0	0	0	0	1	0	1	0	0	1	0	1	0	1	1	1	1	0	0	1	1	0	0	1	1	1	1	0	1	0	0	0	0	0	0	0	0	0	1	0	1	0	0	0	0	1	0	1	1	0	0	1	0	1	0	1	1	1	0	1	1	1];


input=round(rand(1,log2(M)*N)); 
R1=1;%%不要初始相位
 for i=1:3:length(input)-2
    if ((input(i)==1)&(input(i+1)==1)&(input(i+2)==1))
        del_theta=0;
    elseif ((input(i)==0)&(input(i+1)==1)&(input(i+2)==1))
        del_theta=pi/4;
    elseif ((input(i)==0)&(input(i+1)==1)&(input(i+2)==0))
        del_theta=pi/2;
    elseif ((input(i)==0)&(input(i+1)==0)&(input(i+2)==0))
        del_theta=3*pi/4;
    elseif ((input(i)==1)&(input(i+1)==1)&(input(i+2)==0))
        del_theta=-pi/4;
    elseif ((input(i)==1)&(input(i+1)==0)&(input(i+2)==0))
        del_theta=-pi/2;
    elseif ((input(i)==1)&(input(i+1)==0)&(input(i+2)==1))
        del_theta=-3*pi/4;
    elseif ((input(i)==0)&(input(i+1)==0)&(input(i+2)==1))
         del_theta=pi;
     end 
     theta(R1)=del_theta+(R1)*3*pi/8;
     R1=R1+1;
 end
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 x1=exp(j*theta);
 et=x1(1:N);
% plot(x1);
 
%%%%%%%%插3个0=4000%%%%%%%%%%%%%%%%%%%%%
for i=0:length(et)-1
    etx(i*4+1)=real(et(i+1));
    ety(i*4+1)=imag(et(i+1));
    for k=i*4+2:i*4+4
        etx(k)=0;
        ety(k)=0;
    end
end


et_insert_zeros=etx+j*ety;
%gsm/edge pulse shaping  一个符号周期有4个样点的滤波器系数，个数21
%tx=[2.9189e-007	1.888e-005	0.00046755	0.0052743	0.031591	0.11475	0.2828	0.51851	0.75967	0.93598	1	0.93598	0.75967	0.51851	0.2828	0.11475	0.031591	0.0052743	0.00046755	1.888e-005 2.9189e-007];
%%%%%%%%%%%%%%
%测试定时
B=0.3425
R=1;
for t=-T*5/2:Ts:T*5/2   %一个符号间隔内有4个点
h(R)=exp(-1.045*((t-B*T)/T)^2-0.218*((t-B*T)/T)^4);%线性高斯脉冲成形滤波器 加延时了
R=R+1;
end

shaped_et=filter(h,1,[et_insert_zeros]);

figure(2);
plot(shaped_et);
title('pulse shape');
%%%%%%%%%%%%%%%%%%%%%%%%%%%%加高斯噪声%%%%%%%%%%%%%%
I=real(shaped_et);
Q=imag(shaped_et);
gs1=gngauss(sgma);
for i=1:length(I)
  I1(i)=I(i)+gs1;
  Q1(i)=Q(i)+gs1;
end  
shaped_et=I1+j*Q1;


fe=0.001;
Fe=27%频偏110
%Fe=50
thetae0=pi/20%%%%相位pi/32=0.0982,pi/16=0.1963,pi/6=0.5235
%pi/8=0.3927.pi/20=0.1571,pi/25=0.1257

t1=0:Ts:(length(shaped_et)-1)*Ts;
shaped_et=shaped_et.*exp(j*2*pi*Fe*t1);%%%%%%%
shaped_et=shaped_et.*exp(j*thetae0);

figure(3);
plot(shaped_et);
title('gauss + pulse shape+shift');

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%以上为整个加噪声的调制过程

c_opt=[0 0 0 0 0 0 0  -0.012622   0.81641   -5.4860   10.0422   -5.4860    0.81641   -0.012622  0  0  0  0  0  0  0];
%%%21个点的接收滤波器系数
%plot(c_opt);
equalized_st=filter(c_opt,1,[shaped_et]);%zero forcing filter
%%%%%%%%%%定时纠正

k1=1;
I3=real(equalized_st);
Q3=imag(equalized_st);
N3=length(equalized_st);
[I4,Q4,terror1,c_11]=dingshi_modify(I3,Q3,N3);
time=terror1
downsampled_equalizer_output=I4+j*Q4;

%%%%%%%%%%%%55
%%%%%% %%%%%%%%%%%%%%%%%频偏估计加补偿%%%%%%%%%%%%%%%%%%%%%

for i=1:length(downsampled_equalizer_output);
  
    theta2_1(i)=atan2(Q4(i),I4(i));
    theta3_1(i)=theta2_1(i)-(i)*3*pi/8;
  
end
equalized_st1=exp(j*theta3_1);
I1=real(equalized_st1);
Q1=imag(equalized_st1);

[sum,yk,xk,wk]=pin_pian(I1,Q1);
sum1=sum/T;
del_Fe=sum1/(2*pi)
%kalman滤波
%d= [270.833 269.9051 263.5442 259.7248 243.5838 259.4976 247.9547 266.6878 268.9536 260.7016 267.9648 267.4734 259.4608 250.3634 247.0495 267.2073 265.9624 265.8648 245.5509 265.9685];
kal_f = kalmann(del_Fe,  8.8989, 260.7126, 1, 0.01);
del_Fe_kal = kal_f

t=0:T:(length(downsampled_equalizer_output)-1)*T;
downsampled_equalizer_output=downsampled_equalizer_output.*exp(-j*2*pi*del_Fe*t);
%save:
save = downsampled_equalizer_output;%fankui
%%%%%%%%%频偏估计结束%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%相位估计补偿

I5=real(downsampled_equalizer_output);
Q5=imag( downsampled_equalizer_output);
N2=length( downsampled_equalizer_output);
[I6,Q6,theta1]=xiang_wei(I5,Q5,N2);
theta=theta1
downsampled_equalizer_output=I6+j*Q6;
%%%%%%%%%%%%%%%%

downsampled_equalizer_output=downsampled_equalizer_output(6:length(downsampled_equalizer_output));
figure(4);
plot(downsampled_equalizer_output,'.');
%plot(downsampled_equalizer_output);
%title('after receiver with no ISI');

 x=real(downsampled_equalizer_output);
 y=imag(downsampled_equalizer_output);
for i=1:length(downsampled_equalizer_output);
    theta2(i)=atan2(y(i),x(i));
    theta3(i)=theta2(i)-(i)*3*pi/8;
end

x2=exp(j*theta3);
x_judge=real(x2);
y_judge=imag(x2);



for i=1:length(x2)
      phase1(i)=atan2(y_judge(i),x_judge(i));
  end
  %%%%%%%%%%%%%5
figure(5);
plot(exp(j*phase1),'.');
%plot(exp(j*phase1));
  edge_ma_liu1=[];
for i=1:8
    ss(i,:)=[cos(A(i)) sin(A(i))];
end
for i=1:length(phase1)
    r(i,:)=[cos(phase1(i)) sin(phase1(i))];
end

std= 0;  %fankui
for k=1:length(phase1)
    for i=1:8
        cc(i)=dot(r(k,:),ss(i,:));
    end
    cc_max=max(cc);
    if(cc_max==cc(1))
        edge_ma_liu1=[edge_ma_liu1 a4(1,:)];
        std = sqrt( (r(k,1)-ss(1,1))^2 + (r(k,2)-ss(1,2))^2 ) + std;
    elseif(cc_max==cc(2))
        edge_ma_liu1=[edge_ma_liu1 a4(2,:)];
        std = sqrt( (r(k,1)-ss(2,1))^2 + (r(k,2)-ss(2,2))^2 ) + std;
    elseif(cc_max==cc(3))
        edge_ma_liu1=[edge_ma_liu1 a4(3,:)];
        std = sqrt( (r(k,1)-ss(3,1))^2 + (r(k,2)-ss(3,2))^2 ) + std;
    elseif(cc_max==cc(4))
        edge_ma_liu1=[edge_ma_liu1 a4(4,:)];
        std = sqrt( (r(k,1)-ss(4,1))^2 + (r(k,2)-ss(4,1))^2 ) + std;
    elseif(cc_max==cc(5))
        edge_ma_liu1=[edge_ma_liu1 a4(5,:)];
        std = sqrt( (r(k,1)-ss(5,1))^2 + (r(k,2)-ss(5,2))^2 ) + std;
    elseif(cc_max==cc(6))
        edge_ma_liu1=[edge_ma_liu1 a4(6,:)];
        std = sqrt( (r(k,1)-ss(6,1))^2 + (r(k,2)-ss(6,2))^2 ) + std;
    elseif(cc_max==cc(7))
        edge_ma_liu1=[edge_ma_liu1 a4(7,:)];
         std = sqrt( (r(k,1)-ss(7,1))^2 + (r(k,2)-ss(7,2))^2 ) + std;
    elseif(cc_max==cc(8))
        edge_ma_liu1=[edge_ma_liu1 a4(8,:)];
        std = sqrt( (r(k,1)-ss(8,1))^2 + (r(k,2)-ss(8,2))^2 ) + std;
    end
end
std %fankui
m_std =std/length(phase1)%fankui
%---------------------------------------
ps11=0;
for n=1:length(edge_ma_liu1)
    if( input(n)==edge_ma_liu1(n))
        continue;
    else
           ps11=ps11+1;
       end
   end
p11=ps11/length(edge_ma_liu1) 

%----------------------------------------
if (m_std > 0.3)%fankui
    disp('反馈频偏开始');%fankui
    
    pause;%fankui
    p_fankui = fankui(save);%fankui
else%fankui
   %%probility:
ps1=0;
for n=1:length(edge_ma_liu1)
    if( input(n)==edge_ma_liu1(n))
        continue;
    else
           ps1=ps1+1;
       end
   end
p1=ps1/length(edge_ma_liu1) 
%%%%%%%%%%%%%%%
end  %fankui