刘小平课程设计报告.doc

湖南文理学院课程设计报告课程名称： 专业综合课程设计 系 部： 电气与信息工程学院 专业班级： 通信工程07101 学生姓名： 刘小平 （200716020120）指导教师： 侯清莲 完成时间： 2010-6-29 报告成绩： 评阅意见：评阅教师 日期 数字锁相法提取位同步信号设计报告数字锁相法提取位同步信号设计一、设计要求 在任何形式的数字通信系统中，位同步都是必不可少的，无论数字基带传输系统还是数字频带传输系统，无论相干解调还是非相干解调，都必须完成位同步信号的提取，因此，位同步信号的稳定性直接影响到整个数字通信系统的工作性能。根据通信原理所学理论，设计用平方环法从2DPSK信号中提取载波同步信号，并注意相位模糊现象，给出电路结构框图，并完成电路设计、仿真与调试。熟练地运用通信理论，进行数字基带信号、 数字信号频带传输系统、PCM 通信系统和同步系统的设计，并进行通信系统建模。二、设计作用目的（1） 巩固加深载波恢复的认识，提高综合运用通信原理等知识的能力；（2） 培养学生查阅参考文献，独立思考、设计、钻研电子技术相关问题的能力；（3） 通过实际制作安装电子线路，学会单元电路以及整机电路的调试与分析方法；（4） 掌握相关电子线路工程技术规范以及常规电子元器件的性能技术指标；（5） 了解电气图国家标准以及电气制图国家标准，并利用电子CAD正确绘制电路图；（6） 培养严肃认真的工作作风与科学态度，建立严谨的工程技术观念；（7） 培养工程实践能力、创新能力和综合设计能力。三、所用仪器设备计算机 Matlab操作平台四、设计原理4.1 电路分析位同步也称为位定时恢复或码元同步。在任何形式的数字通信系统中，位同步都是必不可少的，无论数字基带传输系统还是数字频带传输系统，无论相干解调还是非相干解调，都必须完成位同步信号的提取，即从接收信号中设法恢复出与发端频率相同的码元时钟信号，保证解调时在最佳时刻进行抽样判决，以消除噪声干扰所导致的解调接收信号的失真，使接收端能以较低的错误概率恢复出被传输的数字信息。因此，位同步信号的稳定性直接影响到整个数字通信系统的工作性能。位同步的实现方法分为外同步法和自同步法两类。由于目前的数字通信系统广泛采用自同步法来实现位同步，故在此仅对位同步中的自同步法进行介绍。采用自同步法实现位同步首先会涉及两个问题：（1）如果数字基带信号中确实含有位同步信息，即信号功率谱中含有位同步离散谱，就可以直接用基本锁相环提取出位同步信号，供抽样判决使用；（2）如果数字基带信号功率谱中并不含有位定时离散谱，怎样才能获得位同步信号。数字基带信号本身是否含有位同步信息与其码型有密切关系。应强调的是，无论数字基带信号的码型如何，数字已调波本身一般不含有位同步信息，因为已调波的载波频率通常要比基带码元速率高得多，位同步频率分量不会落在数字已调波频带之内，通常都是从判决前的基带解调信号中提取位同步信息。二进制基带信号中的位同步离散谱分量是否存在，取决于二进制基带矩形脉冲信号的占空比。若单极性二进制矩形脉冲信号的码元周期为Ts，脉冲宽度为，则NRZ码的 Ts，则NRZ码除直流分量外不存在离散谱分量，即没有位同步离散谱分量1/Ts；RZ码的满足0Ts，且通常占空比为50，此时的RZ码含有n为奇数的n/ Ts离散谱分量，无n为偶数的离散谱分量，这就是说，RZ码含有位同步离散谱分量。显然，为了能从解调后的基带信号中获取位同步信息，可以采取两种措施：（1）如原始数字基带码为NRZ码，若传输信道带宽允许，可将NRZ码变换为RZ码后进行解调；（2）如调制时基带码采用NRZ码，就必须在接收端对解调出的基带信号进行码变换，即将NRZ码变换成RZ码，码变换过程实质上是信号的非线性变换过程，最后再用锁相环（通常为数字锁相环）提取出位同步信号离散谱分量。将NRZ码变为RZ码的最简单的办法是对解调出的基带NRZ码进行微分、整流，即可得到归零窄脉冲码序列。x图1 数字锁相环组成原理框图下面简单介绍一下数字锁相环的组成原理。数字锁相环的主要特点是鉴相信号为数字信号，鉴相输出也是数字信号，即环路误差电压是量化的，没有模拟环路滤波器。由于数字锁相环的输入是经过微分和全波整流后的信号，故这种数字锁相环也称为微分整流型数字锁相环，其原理框图如图1所示。该电路由码型变换器、鉴相器、控制调节器组成，各部分的作用如下： 码型变换器完成解调出的基带NRZ码到RZ码的变换，使鉴相输入信号X含有位同步离散谱分量。 鉴相器用于检测信号X与输出位同步信号（分频输出D）相位间的超前、滞后关系，并以量化形式提供表示实时相位误差的超前脉冲F和滞后脉冲G，供控制调节器使用。当分频输出位同步信号D相位超前与信号X时，鉴相器输出超前脉冲F（低电平有效）；反之，则输出滞后脉冲G（高电平有效），二者均为窄脉冲。 控制调节器的作用是根据鉴相器输出的误差指示脉冲，在信号D与信号X没有达到同频与同相时调节信号D的相位。高稳定晶振源输出180相位差、重复频率为nf0的A、B两路窄脉冲序列作为控制调节器的输入，经n分频后输出重复频率为f0的被调位同步信号D，它与信号X在鉴相器中比相。因超前脉冲F低电平有效并作用于扣除门（与门），平时扣除门总是让脉冲序列A通过，故扣除门为常开门，又因滞后脉冲G高电平有效并作用于附加门（与门），平时附加门总是对序列B关闭的，故附加门为常闭门。当信号D的相位超前与信号X的相位时，鉴相器输出窄的低电平超前脉冲F，扣除门（与门）将从脉冲序列A中扣除一个窄脉冲，则n分频器输出信号D的相位就推迟了Ts /n（相移360/n），信号D的瞬时频率也被调低；当信号D的相位滞后于信号X的相位时，鉴相器输出窄的高电平滞后脉冲G，附加门（与门）此时打开让脉冲序列B（与脉冲序列A保持180固定相差）中的一个脉冲通过，经或门插进来自扣除门输出的脉冲序列A中，则分频器输入多插入的这个脉冲使n分频器输出信号的D相位提前了Ts /n（相移360/n），信号D的瞬时频率则被提高。由此可见，环路对信号D相位和频率的控制调节是通过对n分频器输入脉冲序列步进式加、减脉冲实现的，经环路的这种反复调节，最终可达到相位锁定，从而提取出位同步信号。4.2 性能指标位同步系统的性能通常是用相位误差、建立时间、保持时间等指标来衡量。数字锁相法位同步系统的性能如下。a) 相位误差e数字锁相法提取位同步信号时，相位误差主要是由于位同步脉冲的相位在跳变地调整所引起的。因为每调整一步，相位改变2/n（n是分频器地分频次数），故最大的相位误差为2/n。用这个最大的相位误差来表示，可得 （1）上面已经求得数字锁相法位同步的相位误差有时不用相位差而用时间差Te来表示相位误差。因每码元的周期为T，故得 （2）b) 同步建立时间ts同步建立时间即为失去同步后重建同步所需的最长时间。为了求这个最长时间，令位同步脉冲的相位与输入信号码元的相位相差T/2秒，而锁相环每调整一步仅能移T/n秒，故所需最大的调整次数为 （3）接收随机数字信号时，可近似认为两相邻码元中出现01、10、11、00的概率相等，其中，有过零点的情况占一半。由于数字锁相法中是从数据过零点中提取作比相用的标准脉冲的，因此平均来说，每2T秒可调整一次相位，故同步建立时间为 Ts=2TN=nT（秒） （4）c) 同步保持时间tc当同步建立后，一旦输入信号中断，由于收发双方的固有位定时重复频率之间总存在频差F，收端同步信号的相位就会逐渐发生漂移，时间越长，相位漂移量越大，直至漂移量达到某一准许的最大值，就算失步了。设收发两端固有的码元周期分别为T1=1/F和T2=1/F，则 （5）式中的F0为收发两端固有码元重复频率的几何平均值，且有 （6）由式（5）可得 （7）再由式（6），上式可写为 （8）式（8）说明了当有频差F存在时，每经过T0时间，收发两端就会产生的时间漂移。反过来，若规定两端容许的最大时间漂移为T0/K秒（K为一常数），需要经过多少时间才会达到此值呢？这样求出的时间就是同步保持时间tc。代入式（8）后，得解得 （9）若同步保持时间tc的指标给定，也可由上式求出收发两端振荡器频率稳定度的要求为此频率误差是由收发两端振荡器造成的。若两振荡器的频率稳定度相同，则要求每个振荡器的频率稳定度不能低于 （10）图2 位同步电路原理图只能从码速率为15.625KHz、10KHz、8KHz、4KHz（通过拨码开关SW501选择）的NRZ码中提取出位同步信号。以码速率为15.625KHz的NRZ码为例，将SW501的第一位拨上后，数字锁相环的本振频率就被设置在15.625KHz。在图2中，单片机U508（89C2051）将输入的NRZ码与数字锁相环本振输出的信号的相位进行鉴相（比较两个信号的上升沿），用将相位差进行量化后得到的数值对数字锁相环本振输出的相位进行调整，最后得到正确的位同步信号。4.3、位同步相关知识的简单介绍数字通信中，除了有载波同步的问题外，还有位同步的问题。因为信息是一串相继的信号码元的序列，解调时常需知道每个码元的起止时刻。因此，接收端必须产生一个用作抽样判决的定时脉冲序列，它和接收码元的终止时刻应对齐。我们把在接收端产生与接收码元的重复频率和相位一致的定时脉冲序列的过程称为码元同步或位同步，而称这个定时脉冲序列为码元同步脉冲或位同步脉冲。要使数字通信设备正常工作，离不开正确的位同步信号。如果位同步脉冲发生严重抖动或缺位，则使数字通信产生误码；严重时使通信造成中断。影响位同步恢复的主要原因：输入位同步电路的信号质量；信号的编码方式码元中存在长连“0”或长连“1”。位同步的主要技术指标有静态相差；相位抖动；同步建立时间和同步保持时间。数字通信中位同步恢复的方法主要有两种，一种是发端专门发送导频信号，而另一种是直接从数字信号中提取位同步信号。而直接从数字信号中提取位同步信号也有不止一种方法：滤波法，锁相法两种方法。本课程设计采用的就是用数字锁相环提取位同步信号的方法，这种方法又称为数字锁相。数字锁电路原理图PDLFVCOUiUdUcUo数字锁相原理方框图其工作过程如下：（1）压控振荡器的输出Uo经过采集并分频；（2）输出和基准信号同时输入鉴相器；（3）鉴相器通过比较上述两个信号的频率差，然后输出一个直流脉冲电压Ud；（4）Ud进入到滤波器里面，滤除高频成分后得到信息Ue；（5）Ue进入到压控震荡器VCO里面，控制频率随输入电压线性地变化；（6）这样经过一个很短的时间，VCO的输出就会稳定于某一期望五、具体设计用于提取位同步信号的数字环有超前滞后型数字环和触发器型数字环，系统中的位同步提取模块用的是触发器型数字环，它具有捕捉时间短、抗噪能力强等特点。位同步模块原理框图如图3所示。图3 位同步器方框图位同步模块有以下测试点及输入输出点： S-IN 基带信号输入点/测试点（2个） BS-OUT位同步信号输出点/测试点（3个） 图3中各单元与电路板上元器件的对应关系如下： 晶振 CRY3：晶体；U39：7404 控制器U48：或门7432；U41：计数器74190 鉴相器U40：D触发器7474 量化器U45：可编程计数器8254 数字环路滤波器由软件完成 数控振荡 U46、U45：8254 脉冲展宽器U47：单稳态触发器74123 位同步器由控制器、数字锁相环及脉冲展宽器组成，数字锁相环包括数字鉴相器、量化器、数字环路滤波器、数控振荡器等单元。 下面介绍位同步器的工作原理。 数字锁相环是一个单片机系统，主要器件是单片机89C51及可编程计数器8254。环路中使用了两片8254，共六个计数器，分别表示为8254A0、8254A1、8254A2、8254B0、8254B1、8254B2。它们分别工作在M0、M1、M2三种工作模式。M0为计数中断方式，M1为单稳方式，M2为分频方式。除地址线、数据线外，每个8254芯片还有时钟输入端C、门控信号输入端G和输出端O。数字鉴相器电原理图及波形图如图4（a）、图4（b）所示。输出信号宽度正比于信号ui及uo上升沿之间的相位差，最大值为ui的码元宽度。称此鉴相器为触发器型鉴相器，称包含有触发器型鉴相器的数字环路为触发器型数字锁相环。图4数字鉴相器量化器把相位误差变为多进制数字信号，它由工作于M0方式、计数常数为N0的8254 B2完成（N0为量化级数，此处N0=52）。ud作为8254B2的门控信号，ud为高电平时8254B2进行减计数，ud为低电平时禁止计数，计数结束后从8254B2读得的数字为Nd= N0-Nd 式中Nd为ud脉冲宽度的量化值（下面用量化值表示脉冲宽度和时间间隔），N0Nd，读数结束后再给8254B2写入计数常数N0。读数时刻由8254A2控制，它工作在M1模式，计数常数为N0，ui作为门控信号。一个ui脉冲使8254A2产生一个宽度为N0的负脉冲，倒相后变为正脉冲送到89C51的端，而89C51的外中断1被设置为负跳变中断申请方式。由于8254A2产生的脉冲宽度不小于ud脉冲宽度且它们的前沿处于同一时刻，所以可以确保中断申请后对8254B2读数时它已停止计数。 数字环路滤波器由软件完成。可采用许多种软件算法，一种简单有效的方法是对一组N0作平均处理。设无噪声时环路锁定后ui与uo的相位差为N0/2，则在噪声的作用下，锁定时的相位误差可能大于N0/2也可能小于N0/2。这两种情况出现的概率相同，所以平均处理可以减小噪声的影响，m个Nd值的平均值为 （11）数字滤波器的输出为 Nc = No / 2 + Nd (12）数控振荡器由四个8254计数器及一些门电路构成，其原理框图如图6所示，图中已注明了各个计数器的工作方式和计数常数。 以下分析环路的锁定状态及捕捉过程，此时不考虑噪声的影响。图6 数控振荡器 环路开始工作时，软件使8254B0和8254B1输出高电平，从而使8254A1处于计数工作状态、8254B1处于停止计数状态，G6处于开启状态，8254A1输出一个周期为N0的周期信号。若环路处于锁定状态，则Nd=N0/2，由式（11）及式（12）得Nd=N0/2。此时89c51的P1.4口不输出触发脉冲，8254A0输出端仍保持初始化时的高电平，从而使8254B0的门控端G保持低电平、输出端O保持高电平。这样可保持8254A1、8254B1的工作状态不变、环路仍处于锁定状态。若环路失锁，则NdN0/2，NdN0/2，P1.4口输出一个正脉冲u2，在u2作用下，8254A0输出一个宽度为N0的负脉冲，倒相后变为正脉冲u3送给与门G2。G2的另一个输入信号u1来自8254A1。在G1输出的宽度为N 0的正脉冲持续时间内，8254A1一定有（也只有）一个负脉冲信号输入，此负脉冲经G4倒相后与G1输出的正脉冲相与后给8254B0的G端送一个触发信号u4。在u4的作用下，8254B0输出一个宽度为N0-2的负脉冲。在这段时间内，8254A1停止计数工作，8254B1进行减计数且在此时间内的最后一个时钟周期输出一个负脉冲。8254B0输出的负脉冲的后沿重新启动8254A1，使它重新作N0分频。设m=1，上述过程的有关波形如图7所示，图中uO为环路锁定状态下数控振荡器的输出信号。由图7可见，不管失锁时相位误差多少（不会大于N0），只要对数控振荡器作一次调整，就可使环路进入锁定状态，从而实现快速捕捉。 程序流程如图8所示，输入信号ui使IE1置“1”，且使8254B2计数，对IE1进行位操作时又使之置“0”。由于量化误差，故当Nd为N02，N0/21或N0/21时，环路皆处于锁定状态，不对数控振荡器进行调整。程序中令m=16，进行16次鉴相后做一次平均运算，若发现环路失锁，则对数控振荡器进行一次调整。 控制器的作用是保证每次对8254B2进行读操作之前鉴相器只输出一个正脉冲，它由或门7432（U5:B）及16分频器74190（U13）组成。图7 捕获过程波形当数字环输入信号的码速率与数控振荡器的固有频率完全相同时，环路锁定后输入信号与反馈信号（即位同步信号）的相位关系是固定的且符合抽样判决器的要求（当然开环时它们的相位误差也是固定的，但不符合抽样判决器的要求）。输入信号码速率决定于发送端的时钟频率，数控振荡器固有频率决定于位同步器的时钟频率和数控振荡器固有分频比。由于时钟信号频率稳定度是有限的，故这两个时钟信号的频率不可能完全相同，因此锁相环输入信号码速率与数控振荡器固有频率不可能完全相等（即环路固有频差不为0）。数字环位同步器是一个离散同步器，只有当输入信号的电平发生跳变时才可能对输入信号的相位和反馈信号的相位进行比较从而调整反馈信号的相位，在两次相位调整的时间间隔内，反馈信号的相位相对于输入信号的相位是变化的，即数字环位同步器提取的位同步信号的相位是抖动的，即使输入信号无噪声也是如此。图8 锁相环程序流程 显然，收发时钟频率稳定度越高，数字环的固有频差就越小，提取的位同步信号的相位抖动范围越小。反之，对同步信号的相位抖动要求越严格，则收发时钟的频率稳定度也应越高。 位同步信号抖动范围还与数字位同步器输入信号的连“1”或“0”个数有关，连“1”或“0”个数越多，两次相位调整之间的时间间隔越长，位同步信号的相位抖动越大。 对于NRZ码来说，允许其连“1”、连“0”的个数决定于数字环的同步保持时间tc。tc与收发时钟频率稳定度e、码速率RB、允许的同步误差最大值的关系为：tC =/（2RB）tC的定义是：位同步器输入信号断开后，收发位同步信号相位误差不超过的时间。用模拟环位同步器或模数环位同步器提取的位同步信号的相位抖动与固有频差无关，但随信息码连“1”、连“0”的个数增多而增大。六、仿真结果 输入环路增益为30 仿真时间为5秒 相平面图 输入频率和VCO频率图 输入相位和VCO相位图 频率差图 相位差图 输入环路增益为40 仿真时间为5秒相平面图 输入频率和VCO频率图 相位差图 频率差图 输入频率和VCO频率图 输入相位和VCO相位图 七、心得体会锁相提取位同步信号的设计与仿真，在Matlab的操作环境下实现。该设计涉及到对锁相环的原理、参数选择及应用等各方面的知识，还要求我们对Matlab仿真环境的熟悉和操作能力。锁相环具有载波跟踪特性、调制跟踪特性和低门限特性等优良特性，已经成为电子技术领域中一种相当有效的技术手段，获得了了越来越广泛的应用。在电子技术发展飞速的今天，掌握锁相环的原理及其应用，不管是在学习上还是在以后的工作上，都会产生巨大的作用。学校在强调提高学生的理论水平的时候，也加强了对学生实际动手能力的要求，把学生对理论与实际相结合的操作能力作为了对学生的一项重要考核。由于我们在以前已经有做过Matlab的相关仿真实验和设计，对Matlab的操作环境还是比较熟悉的。而我们在设计的过程中碰到的主要困难是，在给定的锁相环条件下，如何选择合适的参数，以至于在Matlab仿真的条件下，可以得到理想中的结果。锁相技术是这学期开的课，学习的时间不长，对锁相环的了解页仅是基于概念的表面理解。要完成该课程设计，我们就必须深入地理解锁相环的组成、原理及其应用。因此要做的准备工作还是不少。在一个多星期的摸索中，设计完成得不是很好。尽管如此，我还是从中学到了很多东西。通过这个设计，我进一步了解了锁相的相关技术及其应用，也进一步加深了对Matlab仿真的了解，在一定程度上提高了自己的逻辑思维能力和解决问题的能力，有助与加深自己对所学知识的了解和使用，增强了自己的动手能力，为我今后的学习和工作积累了一定的动手实践经验，让我从中受益匪浅。八、参考文献1 郭梯云主编，移动通信，西安电子科技大学，20072 樊昌信主编，通信原理， 国防工业出版社，2007 3 高如云主编，通信电子线路，西安电子科技大学出版社，2007 4吕广平主编，集成电路应用500例，人民邮电出版社，1983九、附录程序代码：% 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; % time if t20 ydd = 4*exp(-t/2)-3*yd*abs(y)-9*y;% de for t=20 end w1b=ydd+w2b;% first integrator - step 1 w2b=ydd+w1b;% first integrator - step 2 yd=w1b/(2*fs); % first integrator output w1c=yd+w2c;% second integrator - step 1 w2c=yd+w1c;% second integrator - step 2 y=w1c/(2*fs); % second integrator output ydv(1,i) = yd;% build dy/dt vector yv(1,i) = y;% build y(t) vector end % 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 = 1 plot(t,fin,k,t,fvco,k) title(Input Frequency and VCO Freqeuncy) xlabel(Time - Seconds);ylabel(Frequency - Hertz);pause elseif k =2 pvco=phin-phierror;plot(t,phin,t,pvco) title(Input Phase and VCO Phase) xlabel(Time - Seconds);ylabel(Phase - Radians);pause elseif k = 3 plot(t,freqerror);title(Frequency Error) xlabel(Time - Seconds);ylabel(Frequency Error - Hertz);pause elseif k = 4 plot(t,phierror);title(Phase Error) xlabel(Time - Seconds);ylabel(Phase Error - Radians);pause elseif k = 5 ppplot elseif k = 6 subplot(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);grid title(Input Frequency and VCO Freqeuncy) xlabel(Time - Seconds);ylabel(Frequency - Hertz);subplot(111) elseif k = 7 kk = 1; end end % End of script file.% File: pllpre.m%clear all % be safedisp( ) % insert blank line fdel = 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;else zeta = 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*nptstsettle = 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 constantsphierror = zeros(1,npts); % initialize vectorfvco=zeros(1,npts);% initialize vector% beginning of simulation loopfor i=1:npts s1=phin(i) - phivco; % phase error s2=sin(s1); % sinusoidal phase detector s3=G*s2; s4=a1*s3; s4a=s4-a2*s5; % loop filter integrator input w1b=s4a+w2b;% filter integrator (step 1) w2b=s4a+w1b;% filter integrator (step 2) s5=w1b/twofs;% generate fiter output s6=s3+s5; % VCO integrator inpu