tran_rec_qdpsk.m

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

%**************           程序说明            ***************************%
%%%%%%%%%%%%%%%% File: tran_rec_dqpsk.m%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%% date: 2007-04-27  author: haige shengzhou corp. %%%%%%%%%%

%%% 目的：进行QDPSK调制解调性能分析 %%%%%%%%%%%%%%%%%%%%%
%%%%载波同步采用costas环
%%%%%%%%%%%%%%%%%%%%%%变量名定义%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

clear all;

nd = 2000;   % Number of symbols that simulates in each loop
 Origin_code_I=rand(1,nd)>0.5;  % 原始I路信号；
 Origin_code_Q=rand(1,nd)>0.5;  % 原始Q路信号；

 [Transimit_code_I,Transimit_code_Q]=qdpskmod(Origin_code_I,Origin_code_Q);%%原始信号编码
ins_value = 8;
for i=1:nd*ins_value
    if rem(i,ins_value) ==1
        Ibit_flowzero(i) =Transimit_code_I(fix((i-1)/ins_value)+1);
    else
        Ibit_flowzero(i) =0;        
    end
end
for i=1:nd*ins_value
    if rem(i,ins_value) ==1
        Qbit_flowzero(i) =Transimit_code_Q(fix((i-1)/ins_value)+1);
    else
        Qbit_flowzero(i) =0;        
    end
end

%升滚降平方根滤波器的参数设置;
alpha=1 ;  %滚降因子。
sr=2.4e3 ; %符号率，符号率为2.4KHz。
%%－－－ins_value倍的采样
fs=sr*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));

 Simulation_Length=length(Transimit_code_I);

  Signal_Source=Transimit_code_I+ j*Transimit_code_Q;

% Signal_Source=Isignal+ j*Qsignal;
%     spow=sum(Isignal.*Isignal+Qsignal.*Qsignal)/nd;  % sum: built in function
% 	attn=0.5*spow*sr/br*10.^(-ebn0/10);
% 	attn=sqrt(attn);  % sqrt: built in function
%     
%     [Transimit_code_I3,Transimit_code_Q3]= comb(Transimit_code_I2,Transimit_code_Q2,attn);% add white gaussian noise
% 	[Transimit_code_I4,Transimit_code_Q4]= compconv(Transimit_code_I3,Transimit_code_Q3,xh2);
%     
%      syncpoint=irfn*IPOINT+1;
%     Transimit_code_I5=Transimit_code_I4(syncpoint:IPOINT:length(Transimit_code_I3));
%     Transimit_code_Q5=Transimit_code_Q4(syncpoint:IPOINT:length(Transimit_code_Q3));
offset=0; % 
ideal_samplefre =ins_value*sr;
samplefre= ideal_samplefre+offset;
%－－－－－－－－生成有载波调制的数据－－－－－－－－－－－－－－－－－－－－－－
fcarrier = 2e5;
foffset = 0;
Freq_Sample=6e5;
ophase= pi/3;
time=[0:length(Transimit_code_I)-1]/Freq_Sample;
Carrier=exp(j*(fcarrier/Freq_Sample*(1:Simulation_Length)+ophase));
% time=[0:length(Isignal)-1]/Freq_Sample;
% Carrier=exp(j*(fcarrier/Freq_Sample*(1:Simulation_Length)+ophase));

Signal_Channel=Signal_Source.*Carrier;
% Signal_Channel=Transimit_code_I.*cos((fcarrier/Freq_Sample*(1:Simulation_Length)+ophase))-Transimit_code_I.*sin((fcarrier/Freq_Sample*(1:Simulation_Length)+ophase));

Signal_PLL=zeros(Simulation_Length,1);
NCO_Phase = ones(Simulation_Length,1);
Discriminator_Out=zeros(Simulation_Length,1);
Freq_Control=zeros(Simulation_Length,1);
PLL_Phase_Part=zeros(Simulation_Length,1);
PLL_Freq_Part=zeros(Simulation_Length,1);
%环路处理
Bd=1000;                              % 环路噪声带宽
damp=0.707;                          % 阻尼系数
td=1/10e3;                           % 相关积分时间
Kd=1;                                % 环路增益
Wd=2*Bd/(damp+1/(4*damp));           % 自然角频率 
C1=8*damp*Wd*td/(Kd*(4+4*damp*Wd*td+(Wd*td)^2));  % 环路滤波器的系数
C2=4*(Wd*td)^2/(Kd*(4+4*damp*Wd*td+(Wd*td)^2));
for i=2:Simulation_Length

Signal_PLL(i)=Signal_Channel(i)*exp(-j*mod(NCO_Phase(i-1),2*pi));
I_PLL(i)=real(Signal_PLL(i));
Q_PLL(i)=imag(Signal_PLL(i));
% I_PLL(i)=Signal_Channel(i)*cos(mod(NCO_Phase(i-1),2*pi));
% Q_PLL(i)=-1*Signal_Channel(i)*sin(mod(NCO_Phase(i-1),2*pi));

% if PLL_Mode == 1
% Discriminator_Out(i)=atan2(Q_PLL(i),I_PLL(i));
% else if PLL_Mode == 2
% Discriminator_Out(i)=sign(I_PLL(i))*Q_PLL(i)/abs(Signal_PLL(i));
% else
Discriminator_Out(i)=(sign(I_PLL(i))*Q_PLL(i)-sign(Q_PLL(i))*I_PLL(i))...
/(sqrt(2)*sqrt(I_PLL(i)^2+Q_PLL(i)^2));
PLL_Phase_Part(i)=Discriminator_Out(i)*C1;
Freq_Control(i)=PLL_Phase_Part(i)+PLL_Freq_Part(i-1);
PLL_Freq_Part(i)=Discriminator_Out(i)*C2+PLL_Freq_Part(i-1);
NCO_Phase(i)=NCO_Phase(i-1)+Freq_Control(i);
end

figure
subplot(1,2,1)
plot(PLL_Freq_Part(2:Simulation_Length)*Freq_Sample);
grid on;
title('锁相环频率响应曲线');
% axis([1 Simulation_Length -100 100]);
subplot(1,2,2)
plot(PLL_Phase_Part(2:Simulation_Length)*180/pi);
title('锁相环相位响应曲线');
% axis([1 Simulation_Length -2 2]);
grid on;
figure(2)
plot(I_PLL);
figure(3)
plot(Q_PLL);
figure(4)
subplot(2,1,1)
plot(I_PLL(1500:1550));
subplot(2,1,2)
plot(Transimit_code_I(1500:1550));
figure(5)
subplot(2,1,1)
plot(I_PLL(1500:1550));
subplot(2,1,2)
plot(-Transimit_code_I(1500:1550));
       
%-----------------------------接收机部分－－－－－－－－－－－－－
% fc = 21.4e6;
% LP_fircof = fir1(128,0.125);
%     Rec_Ibasesignal =tra_IFsignal.*cos(2*pi*fc.*time+ophase);
%     Rec_Qbasesignal=-1*tra_IFsignal.*sin(2*pi*fc.*time+ophase);
%     %－－－－－－－－－－－－－－－匹配滤波－－－－－－－－
%     Rec_MFIbase =filter(psf,1,Rec_Ibasesignal);%产生(length(LP_fircof)-1)/2位延迟
%         Rec_MFQbase =filter(psf,1,Rec_Qbasesignal);         %－－－－－－－－－－－－－－－低通滤波－－－－－－－－
%         Rec_LPIbase =filter(LP_fircof,1,Rec_MFIbase);%产生(length(LP_fircof)-1)/2位延迟
%         Rec_LPQbase=filter(LP_fircof,1,Rec_MFQbase); 
%         FIRfilter_delay =(length(LP_fircof)-1)/2+(length(psf)-1)/2;
%         Rec_Idec= Rec_LPIbase(FIRfilter_delay+1:ins_value:length(Rec_LPIbase));
%         Rec_Qdec= Rec_LPQbase(FIRfilter_delay+1:ins_value:length(Rec_LPQbase));
%     
%     Transimit_code_I6=Rec_Idec>=0;
%     Transimit_code_Q6=Rec_Qdec>=0;
%     [receive_decode_I,receive_decode_Q]=dpskdemod(Transimit_code_I6,Transimit_code_Q6);%译码

%     
%     figure(1)
%     subplot(2,1,1);
%     plot(Origin_code_I);
%      subplot(2,1,2);
%     plot(receive_decode_I);

%      figure(2)
%     subplot(2,1,1);
%     plot(Transimit_code_I);
%      subplot(2,1,2);
%     plot(Transimit_code_I5);
%        
%       figure(3)
%     subplot(2,1,1);
%     plot(Transimit_code_Q);
%      subplot(2,1,2);
%     plot(Transimit_code_Q5);
%************************** Bit Error Rate (BER) ****************************
% Transimit_code_I6=2*Transimit_code_I6-1;
% Transimit_code_Q6=2*Transimit_code_Q6-1;
%     noe1=sum(abs(Transimit_code_I6-Transimit_code_I))+sum(abs(Transimit_code_Q6-Transimit_code_Q));
% 
%     noe2=sum(abs(receive_decode_I-Origin_code_I))+sum(abs(receive_decode_Q-Origin_code_Q));  % sum: built in function
% 	nod2=length(Origin_code_I)+length(Origin_code_Q);  % length: built in function
% 	noe=noe+noe2;
% 	nod=nod+nod2;
% 
% end % for iii=1:nloop    
% 
% %********************** Output result ***************************
% 
% ber = noe/nod;