bpskcostas.m

实现QPSK和BPSK的调制解调程序

%**************           程序说明            ***************************%
%%%%%%%%%%%%%%%% File: bpskcostas.m%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%% date: 2007-06-14  author: haige shengzhou corp. %%%%%%%%%%

%%% 目的：进行BPSK调制性能分析 %%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%变量名定义%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%程序主体%%%%%%%%%%%%%%%%%%%%%%%%%
clear all;
K=1;%环路增益
phase=0:359;
cosdata=cos(phase/180*pi);%%三角函数对应角度值
sindata=sin(phase/180*pi);
realdatafifo=zeros(1,32);
imagdatafifo=zeros(1,32);
errorfifo=zeros(1,32);
out_phaseerror=0;
nd = 500;   % Number of symbols that simulates in each loop
for i=1:nd
    if rand>0.5
        trans_code(i)=1;
      Origin_code_I(i)=cos(pi);
      Origin_code_Q(i)=sin(pi);
  else
      trans_code(i)=0;
      Origin_code_I(i)=cos(0);
      Origin_code_Q(i)=sin(0);
  end
end
ins_value = 16;
for i=1:nd*ins_value
    i
    if rem(i,ins_value) ==1
        Ibit_flowzero(i) =Origin_code_I(fix((i-1)/ins_value)+1);
    else
        Ibit_flowzero(i) =0;        
    end
end
for i=1:nd*ins_value
    i
    if rem(i,ins_value) ==1
        Qbit_flowzero(i) =Origin_code_Q(fix((i-1)/ins_value)+1);
    else
        Qbit_flowzero(i) =0;        
    end
end
%升滚降平方根滤波器的参数设置;
alpha=1 ;  %滚降因子。
sr=1e3 ; %符号率，符号率为10KHz。内插1000，数据率为10M;
fs=6e6;
ts=1/fs;
fc=2e6;%载波频率
fset=0;
%%－－－ins_value倍的采样

[psf,den] = rcosine(1*sr,ins_value*sr,'fir/normal',alpha,1);
Ipulse=conv(Ibit_flowzero,psf);
Qpulse=conv(Qbit_flowzero,psf);% 卷积后数据长度为length(Qbit_flowzero)+length(psf)-1
%去处滤波器延迟效应
Isignal =Ipulse((length(psf)-1)/2+1:length(Ipulse));
Qsignal =Qpulse((length(psf)-1)/2+1:length(Qpulse));
%－－－－－－－－生成有载波调制的数据－－－－－－－－－－－－－－－－－－－－－－
ophase= 0.3*pi;
time=[0:length(Isignal)-1]/fs;
tra_IFsignal=Isignal.*cos(2*pi*(fc+fset).*time+ophase)-Qsignal.*sin(2*pi*(fc+fset).*time+ophase);
%%%%-------解调载波同步采用costas环－-------------------------------------------
LPF_coffic=fir1(31,0.2);%%%低通滤波器参数
H_coffic=fir1(31,0.1);%%%环路滤波器参数
outphase=0;%fix(rand*10);%%初始相位随机
origi_phaseerror=0;
phased=fix(fc/fs*360);
cosvalue=cosdata(outphase+1);
sinvalue=sindata(outphase+1);
errorfreqold=0;

    Rec_Ibasesignal =tra_IFsignal.*cos(2*pi*fc.*time+ophase);
    Rec_Qbasesignal=tra_IFsignal.*sin(2*pi*fc.*time+ophase);
    
    phin(1)=rand*pi;
for n=2:length(Isignal)
    n
    
    realdata(n)=tra_IFsignal(n)*cos(phin(n-1));
    imagdata(n)=tra_IFsignal(n)*sin(phin(n-1));
    phaseoffset=0;
    for m=32:-1:2
        realdatafifo(m)=realdatafifo(m-1);
        imagdatafifo(m)=imagdatafifo(m-1);
    end
    realdatafifo(1)=realdata(n);
    imagdatafifo(1)=imagdata(n);%%%采用倒序便于卷积滤波
    
    realerror=sum(realdatafifo.*LPF_coffic);
    realerrorout(n)=realerror;
    imagerror=sum(imagdatafifo.*LPF_coffic);
    imagerrorout(n)=imagerror;
    errorfreqnew=atan2(imagerror,realerror);
    errorf(n)=(errorfreqnew-errorfreqold)/ts/(2*pi);
    errorfreqold=errorfreqnew;
    
    phase_error=realerror*imagerror;
    for m=32:-1:2
        errorfifo(m)=errorfifo(m-1);
    end
    errorfifo(1)=phase_error;
    out_phaseerror=sum(errorfifo.*H_coffic);
    phaseoffsetout(n)=out_phaseerror;
%     while mod(n,16)==0;
%     out_phaseerror=sum(errorfifo.*H_coffic);
%     phaseoffset=fix(K*out_phaseerror*100);%%%%%相位偏移值，在得到后只改变一次
%     break;
%     end;
%     phaseoffsetout(n)=phaseoffset;
%     
%     outerror(n)=out_phaseerror;
%     outphase=outphase+phased+phaseoffset;
%      if outphase>=360
%         outphase=outphase-360;
%     end
%      if outphase<0
%          outphase=outphase+360;
%      end
%      if outphase==0
%          outphase=outphase+1;
%      end
%   
%     cosvalue=cosdata(outphase);
%     sinvalue=sindata(outphase);
%     outcosdata(n)=cosvalue;
%     outsindata(n)=sinvalue;
phin(n)=phin(n-1)+2*pi*(fc)*ts+50*out_phaseerror;

end
recive_code=realerrorout(1:ins_value:length(Isignal));
recive_decode=recive_code<0;
plotcode1=trans_code(400:500);
plotcode2=recive_decode(400:500);
figure(1);
plot(errorf);
figure(2)
subplot(2,1,1)
plot(-10*realerrorout(2000:length(Isignal)));
subplot(2,1,2)
plot(Isignal(2000:length(Isignal)));
figure(3)
subplot(2,1,1)
plot(imagerrorout(2000:length(Isignal)));
subplot(2,1,2)
plot(Qsignal(2000:length(Isignal)));
figure(3)
plot(realerrorout);
figure(4)
plot(phaseoffsetout);
figure(5)
subplot(2,1,1)
plot(plotcode1);
subplot(2,1,2)
plot(plotcode2);

figure(6)
subplot(2,1,1)
plot(plotcode1);
subplot(2,1,2)
plot(-1*plotcode2);