pllsim.m

二阶锁相环Matlab仿真代码

%PLLSim.m
%锁相环仿真
%使用QPSK调制，标准锁相环

clear;                                       %清理内存空间
clc;                                         %清屏
tic                                          %时间统计开始
close;
%%%%%%%%%%仿真参数设置%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%        aa 
DataModulation=1;							%=1:有符号调制，=0：无符号调制，一直发+1
AddNoice=1;									%=1:加入噪声，=0：完全不加噪声
IsQSPK=0;                                   %=0: BSPK, =1:QPSK
CarrierErrOn=1;                             %=1:加入频偏，=0：不加入频偏
PllOn=1;                                    %=1: use Pll, =0:don't use pll

snr_in_dB=10;                                 %Eb/N0 value expressed in dB
SymNum=10000;                                 %Nomber of symbols to be simulated.

%频偏参数
InitFeqErr=-0/64;     							%初始频偏，以符号速率归一化，即频偏是符号速率的df倍
FeqErrAcc=0.0000;									%频差的斜生率，以符号速率归一化，即频偏是符号速率的df倍
MaxFeqErr=120/64;									%最大的频偏值
InitPhase=0.1*pi;								%载波初始相位，单位:rad

%锁相环参数
Pll_fb_a=1;
Pll_b=-0.9554388;
Pll_A=1;
Pll_K=0.1782447;

%%%%%%%%%%约定参数与参数换算%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if IsQSPK==1
	M=4;                                         %调制阶数，including log2(M) bits per symbol.
else
	M=2;
end
Eb=1;                          		  		 %Energy of per bit.
SNR=10^(snr_in_dB/10);                       %Converted to expressed by times(以倍数表示的Eb/N0).

%%%%%%%%%%%仿真初始化%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
	%调制
	if DataModulation==1 					%=1:有符号调制，=0：无符号调制，一直发+1
	    dsource_i=1-2*randint(1,SymNum);    %Date source
	    dsource_q=1-2*randint(1,SymNum);	%Qpsk modulation. Map 0->+1, 1->-1.
    else
    	dsource_i=ones(1,SymNum);
    	dsource_q=ones(1,SymNum);
    end
    if IsQSPK==1
    	SendSig=dsource_i+j.*dsource_q;		%发送的信号（复数）
    else
    	SendSig=dsource_i;
    end
    SymIdx=1:SymNum;
    SendSig=SendSig.*exp(j*10.*SymIdx);		%载波调制，即使是BPSK也要调制，以便作为复数传给AWGN信道
    %信道
   	if AddNoice==1								%=1:加入噪声，=0：完全不加噪声
		RecSig=awgn(SendSig,snr_in_dB);			%
   	else
   		RecSig=SendSig;
   	end
   	RecSig=RecSig.*exp(-j*10.*SymIdx);		%下变频
   	
   	if CarrierErrOn==1								%=1:加入频偏，=0：不加入频偏
   		FeqErr(1)=InitFeqErr;						%初始频差
   		CarrierPhase(1)=InitPhase;					%初始相差
   		for m=1:SymNum
   			RecSig(m)=RecSig(m)*exp(j*CarrierPhase(m));		%加入载波频偏
   			FeqErr(m+1)=FeqErr(m)+FeqErrAcc;            	%频差斜升
   			if abs(FeqErr(m+1))>=MaxFeqErr					%达到最大频偏时翻转,成为锯齿形状
   				FeqErrAcc=-FeqErrAcc;
   			end
   			CarrierPhase(m+1)=CarrierPhase(m)+FeqErr(m+1)*2*pi;	%载波相位累加
   		end
   		subplot(2,2,1);
   		plot(FeqErr(1:SymNum));
   	end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%   	
   	%接收
   	if PllOn==1
	   	FeqAccum=0;   %环路滤波器初始累加值
	   	NCOPhase(1)=0;		%NCO相位初值
	   	for SymIdx=1:SymNum
	   		NCOSig(SymIdx)=RecSig(SymIdx)*exp(j*NCOPhase(SymIdx));
	   	
	   		%相差检测
		   	if IsQSPK==0  %BPSK
		   		PhaseErr=-sign((real(NCOSig(SymIdx))))*imag(NCOSig(SymIdx));
		   		%仿真比较表明sign(cos)*sin比cos*sin 鉴相更快，大约1倍
%		   		PhaseErr=(real(NCOSig(SymIdx)))*imag(NCOSig(SymIdx));
            else  %QPSK
            	PhaseErr=-sign((real(NCOSig(SymIdx))))*imag(NCOSig(SymIdx))+sign(imag(NCOSig(SymIdx)))*(real(NCOSig(SymIdx)));
		   	end
		   	
	%	   	环路滤波器
		   	PhaseErrAK=PhaseErr*Pll_A*Pll_K;
		   	FeqAccum=FeqAccum*Pll_fb_a+(1+Pll_b)*PhaseErrAK;
		   	Fcw=FeqAccum+PhaseErrAK;
	%	   	NCO
			NCOPhase(SymIdx+1)=NCOPhase(SymIdx)+Fcw;
		end
		%	对相位偏差画图
			subplot(2,2,3);
			PhaseDiff=NCOPhase(1:SymNum)+CarrierPhase(1:SymNum);
			ModValue=2*pi/M;   			 %模糊度的值
			PhaseDiff=mod(PhaseDiff,ModValue);
			PhaseIdx=find(PhaseDiff>ModValue/2);
			PhaseDiff(PhaseIdx)=PhaseDiff(PhaseIdx)-ModValue;
			plot(PhaseDiff);
	else
		NCOSig=RecSig;
	end
	subplot(2,2,2);
   	plot(real(NCOSig));
	subplot(2,2,4);
	plot(imag(NCOSig));
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%     	
   	%判决
   	numofbiterror=0;   
   	DetectSig=NCOSig;
   	SynLength=200;				%同步所需符号数，即统计误码时需要扣除的符号数
   	SignDet_i=sign(real(DetectSig));
    numofbiterror=numofbiterror+sum(SignDet_i(SynLength+1:end)~=dsource_i(SynLength+1:end));
    if IsQSPK==1
   	 	SignDet_q=sign(imag(DetectSig));
   		numofbiterror=numofbiterror+sum(SignDet_q(SynLength+1:end)~=dsource_q(SynLength+1:end));
	end   

	ber=numofbiterror/(log2(M)*(SymNum-SynLength))  %Calculate the bit error rate. Has log2(M)*N bits.
	theo_err_ber=Qfunct(sqrt(2*SNR))             %BPSK/QPSK Theory bit error rate.
toc                                          %时间统计结束
%%=========================