锁相环仿真(基于MATLAB)

锁相环仿真1.锁相环的理论分析1．1锁相环的基本组成锁相环路是一种反馈控制电路，简称锁相环（PLL,Phase-Locked Loop）。

锁相环的特点是：利用外部输入的参考信号控制环路内部振荡信号的频率和相位。

因锁相环可以实现输出信号频率对输入信号频率的自动跟踪，所以锁相环通常用于闭环跟踪电路。

锁相环在工作的过程中，当输出信号的频率与输入信号的频率相等时，输出电压与输入电压保持固定的相位差值，即输出电压与输入电压的相位被锁住，这就是锁相环名称的由来。

锁相环通常由鉴相器（PD,Phase Detector）、环路滤波器（LF,Loop Filter）和压控振荡器（VCO,Voltage Controlled Oscillator）三部分组成，锁相环组成的原理框图如图示：锁相环中的鉴相器又称为相位比较器，它的作用是检测输入信号和输出信号的相位差，并将检测出的相位差信号转换成u（t）电压信号输出，该信号经低通滤D（t），对振荡器输出信号的频率实施波器滤波后形成压控振荡器的控制电压uC控制。

1．2锁相环的工作原理1.2.1鉴相器锁相环中的鉴相器（PD）通常由模拟乘法器组成，利用模拟乘法器组成的鉴相器电路如图示：鉴相器的工作原理是：设外界输入的信号电压和压控振荡器输出的信号电压分别为：为压控振荡器在输入控制电压为零或为直流电压时的振荡角频率，称式中的ω为：为电路的固有振荡角频率。

则模拟乘法器的输出电压uD1.2.2 低通滤波器低通滤波器（LF）的将上式中的和频分量滤掉，剩下的差频分量作为压控振荡器的输入控制电压uC （t）。

即uC（t）为：式中的ωi 为输入信号的瞬时振荡角频率，θi（t）和θO（t）分别为输入信号和输出信号的瞬时位相，根据相量的关系可得瞬时频率和瞬时位相的关系为：即则，瞬时相位差θd为对两边求微分，可得频差的关系式为上式等于零，说明锁相环进入相位锁定的状态，此时输出和输入信号的频率和相位保持恒定不变的状态，uc（t）为恒定值。

当上式不等于零时，说明锁相环的相位还未锁定，输入信号和输出信号的频率不等，uc（t）随时间而变。

1.2.3 压控振荡器压控振荡器（VCO）的压控特性如图示该特性说明压控振荡器的振荡频率ωu 以ω为中心，随输入信号电压uc（t）线性地变化，变化的关系如下：上式说明当uc （t）随时间而变时，压控振荡器（VCO）的振荡频率ωu也随时间而变，锁相环进入“频率牵引”，自动跟踪捕捉输入信号的频率，使锁相环进入锁定的状态，并保持ω0=ωi的状态不变。

2．信号流程图锁相环的原理框图如下：其工作过程如下：（1）压控振荡器的输出Uo经过采集并分频；（2）输出和基准信号同时输入鉴相器；（3）鉴相器通过比较上述两个信号的频率差，然后输出一个直流脉冲电压Ud；（4）Ud进入到滤波器里面，滤除高频成分后得到信息Ue；（5）Ue进入到压控震荡器VCO里面，控制频率随输入电压线性地变化；（6）这样经过一个很短的时间，VCO的输出就会稳定于某一期望值。

3．二阶环仿真源程序代码及仿真结果3．1 程序代码：% File: c6_nltvde.mw2b=0; w2c=0; % initialize integratorsyd=0; y=0; % initialize differential equationtfinal = 50; % simulation timefs = 100; % sampling frequencydelt = 1/fs; % sampling periodnpts = 1+fs*tfinal; % number of samples simulatedydv = zeros(1,npts); % vector of dy/dt samplesyv = zeros(1,npts); % vector of y(t) samples%% beginning of simulation loopfor i=1:nptst = (i-1)*delt; % timeif t<20ydd = 4*exp(-t/2)-3*yd*abs(y)-9*y; % de for t<20elseydd = 4*exp(-t/2)-3*yd-9*y; % de for t>=20endw1b=ydd+w2b; % first integrator - step 1w2b=ydd+w1b; % first integrator - step 2yd=w1b/(2*fs); % first integrator outputw1c=yd+w2c; % second integrator - step 1 w2c=yd+w1c; % second integrator - step 2 y=w1c/(2*fs); % second integrator outputydv(1,i) = yd; % build dy/dt vectoryv(1,i) = y; % build y(t) vectorend % end of simulation loopplot(yv,ydv) % plot phase planexlabel('y(t)') % label x axisylabel('dy/dt') % label y zxis% End of script file.% File: pllpost.m%kk = 0;while kk == 0k = menu('Phase Lock Loop Postprocessor',...'Input Frequency and VCO Frequency',...'Input Phase and VCO Phase',...'Frequency Error','Phase Error','Phase Plane Plot',...'Phase Plane and Time Domain Plots','Exit Program');if k == 1plot(t,fin,'k',t,fvco,'k')title('Input Frequency and VCO Freqeuncy')xlabel('Time - Seconds');ylabel('Frequency - Hertz');pause elseif k ==2pvco=phin-phierror;plot(t,phin,t,pvco)title('Input Phase and VCO Phase')xlabel('Time - Seconds');ylabel('Phase - Radians');pause elseif k == 3plot(t,freqerror);title('Frequency Error')xlabel('Time - Seconds');ylabel('Frequency Error - Hertz');pause elseif k == 4plot(t,phierror);title('Phase Error')xlabel('Time - Seconds');ylabel('Phase Error - Radians');pause elseif k == 5ppplotelseif k == 6subplot(211);phierrn = phierror/pi;plot(phierrn,freqerror,'k');grid;title('Phase Plane Plot');xlabel('Phase Error /Pi');ylabel('Frequency Error - Hertz');subplot(212)plot(t,fin,'k',t,fvco,'k');gridtitle('Input Frequency and VCO Freqeuncy')xlabel('Time - Seconds');ylabel('Frequency - Hertz');subplot(111) elseif k == 7kk = 1;end end % End of script file.% File: pllpre.m%clear all % be safedisp(' ') % insert blank linefdel = input('Enter the size of the frequency step in Hertz > ');fn = input('Enter the loop natural frequency in Hertz > ');lambda = input('Enter lambda, the relative pole offset > ');disp(' ')disp('Accept default values:')disp(' zeta = 1/sqrt(2) = 0.707,')disp(' fs = 200*fn, and')disp(' tstop = 1')dtype = input('Enter y for yes or n for no > ','s');if dtype == 'y'zeta = 1/sqrt(2);fs = 200*fn;tstop = 1;elsezeta = input('Enter zeta, the loop damping factor > ');fs = input('Enter the sampling frequency in Hertz > ');tstop = input('Enter tstop, the simulation runtime > ');end %npts = fs*tstop+1; % number of simulation pointst = (0:(npts-1))/fs; % default time vectornsettle = fix(npts/10); % set nsettle time as 0.1*npts tsettle = nsettle/fs; % set tsettle% The next two lines establish the loop input frequency and phase% deviations.fin = [zeros(1,nsettle),fdel*ones(1,npts-nsettle)];phin = [zeros(1,nsettle),2*pi*fdel*t(1:(npts-nsettle))];disp(' ') % insert blank line% end of script file pllpre.m % File: pll2sin.mw2b=0; w2c=0; s5=0; phivco=0; %initializetwopi=2*pi; % define 2*pitwofs=2*fs; % define 2*fsG=2*pi*fn*(zeta+sqrt(zeta*zeta-lambda)); % set loop gaina=2*pi*fn/(zeta+sqrt(zeta*zeta-lambda)); % set filter parametera1=a*(1-lambda); a2 = a*lambda; % define constants phierror = zeros(1,npts); % initialize vectorfvco=zeros(1,npts); % initialize vector % beginning of simulation loopfor i=1:nptss1=phin(i) - phivco; % phase errors2=sin(s1); % sinusoidal phase detectors3=G*s2;s4=a1*s3;s4a=s4-a2*s5; % loop filter integrator inputw1b=s4a+w2b; % filter integrator (step 1)w2b=s4a+w1b; % filter integrator (step 2)s5=w1b/twofs; % generate fiter outputs6=s3+s5; % VCO integrator inputw1c=s6+w2c; % VCO integrator (step 1)w2c=s6+w1c; % VCO integrator (step 2)phivco=w1c/twofs; % generate VCO outputphierror(i)=s1; % build phase error vectorfvco(i)=s6/twopi; % build VCO input vectorend% end of simulation loopfreqerror=fin-fvco; % build frequency error vector % End of script file.function [] = pplane(x,y,nsettle)% Plots the phase plane with phase in the range (-pi,pi)ln = length(x);maxfreq = max(y);minfreq = min(y);close % Old figure discardedaxis([-1 1 1.1*minfreq 1.1*maxfreq]); % Establish scalehold on % Collect info for new figj = nsettle;while j < lni = 1;while x(j) < pi & j < lna(i) = x(j)/pi;b(i) = y(j);j = j+1;i = i+1;endplot(a,b,'k')a = [];b = [];x = x - 2*pi;endhold offtitle('Phase-Plane Plot')xlabel('Phase Error / Pi')ylabel('Frequency Error in Hertz')grid % End of script file.% File: ppplot.m% ppplot.m is the script file for plotting phase plane plots. If the% phase plane is constrained to (-pi,pi) ppplot.m calls pplane.m.kz = 0;while kz == 0k = menu('Phase Plane Options',...'Extended Phase Plane',...'Phase Plane mod(2pi)',...'Exit Phase Plane Menu');if k == 1phierrn = phierrn/pi;plot(phierrn,freqerror,'k')title('Phase Plane Plot')xlabel('Phase Error /Pi')ylabel('Frequency Error - Hertz')gridpauseelseif k == 2pplane(phierrn,freqerror,nsettle+1)pauseelseif k == 3kz = 1;endend % End of script file.3．2仿真结果：G=30时的仿真图形：Accept the tentative values:the first loop frequency is 5Enter y for yes or n for no > yEnter the loop gain >30 输入环路增益为30Enter the sampling frequency in Hertz > 1200Enter tstop, the simulation runtime > 5 仿真时间为5秒相平面图输入频率和VCO频率图输入相位和VCO相位图频率差图相位差图G=40时的仿真图形：Accept the tentative values:the first loop frequency is 5Enter y for yes or n for no > yEnter the loop gain >40 输入环路增益为40Enter the sampling frequency in Hertz > 1200Enter tstop, the simulation runtime > 5 仿真时间为5秒相平面图输入频率和VCO频率图相位差图频率差图输入频率和VCO频率图输入相位和VCO相位图。