czy.m

在Windows 操作系统下

clear all
disp('请输入： ')						% 删除所有以前的数据
disp('提示：输入阶跃信号频率和固有频率应相差不大')
fdel = input('输入阶跃信号频率（Hz）> ');
fn = input('固有频率（Hz）> ');
disp('相对极点偏移： ')
lambda = input('要求小于阻尼因子zeta（默认值：0.707）的平方 > ');

disp(' ')
disp(' 以下为系统默认值:')
disp(' 阻尼因子zeta = 1/sqrt(2) = 0.707,')
disp(' 抽样频率fs = 20*fn')
disp(' 仿真停止时间（s）tstop = 1')
disp(' ')

disp(' 请选择：')
dtype = input('接受默认值请输入：y,自行设置参数请输入：n > ','s');
if dtype == 'y'
   zeta = 1/sqrt(2);
   fs = 20*fn;
   tstop = 1;
else 
disp(' ')
disp('请输入： ')	
zeta = input('锁相环的阻尼因子 > ');
fs = input('抽样频率（Hz） > ');
tstop = input('仿真运行的时间（s） > ');
end


npts = fs*tstop+1;                 % 仿真的点数
t = (0:(npts-1))/fs;            % 时间向量
nsettle = fix(npts/10);         % 设置0.1*npts时间停留
tsettle = nsettle/fs;           % 设置开始时间
%
%下面一行是进入环路的相位
%
fin = [zeros(1,nsettle),fdel*ones(1,npts-nsettle)];
phin = [zeros(1,nsettle),2*pi*fdel*t(1:(npts-nsettle))];

disp(' ')                      

w2b=0; w2c=0; s5=0; phivco=0;				%初始化    
twopi=2*pi;									% 定义2*pi,pi的倍
twofs=2*fs;									% 定义2*fs,抽样频率的倍
G=2*pi*fn*(zeta+sqrt(zeta*zeta-lambda));	% 设置环路鉴相器的增益
a=2*pi*fn/(zeta+sqrt(zeta*zeta-lambda));	% 设置环路滤波器的参数
a1=a*(1-lambda); a2 = a*lambda;				% 定义常量,环路滤波器的指标
phierror = zeros(1,npts);	  				% 初始化相位差矢量,这里是全0
fvco=zeros(1,npts);							% 

%
% 开始仿真锁相环
%

for i=1:npts
        s1=phin(i) - phivco;				% 计算相位差
        s2=sin(s1);	                        % 由相位差产生正弦特性的误差电压
        s3=G*s2;
        s4=a1*s3;
        s4a=s4-a2*s5;                   	% 环路滤波器的输入
        w1b=s4a+w2b;						% 滤波器的第1步
        w2b=s4a+w1b;						% 滤波器的第2步					
        s5=w1b/twofs;						% 滤波器的输出
        s6=s3+s5;                       	% VCO 的输入
        w1c=s6+w2c;							% VCO 第1步综合
        w2c=s6+w1c;							% VCO 第2步综合
        phivco=w1c/twofs;					% 产生 VCO 输出
        phierror(i)=s1;						% 存储相位误差
        fvco(i)=s6/twopi;					% VCO输出矢量
end
%
% 仿真结束
%
freqerror=fin-fvco;								% 计算频率差
% 以下画图
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
figure
subplot(2,2,1)
plot([1:npts],fin);
title('锁相环输入频率曲线');

subplot(2,2,2)
plot(fvco);
title('锁相环输出频率曲线');

subplot(2,2,3)
plot(freqerror);
title('锁相环频率差曲线');

subplot(2,2,4)
plot(phierror);
title('锁相环相位差曲线');

figure
subplot(2,1,1)
%plot(t,sin(fvco*2*pi+phierror));
plot(t,sin((fvco*2*pi+phierror).*t));
title('锁相环的输出波形');

subplot(2,1,2)
plot(t,sin(fin*2*pi.*t));
title('进入锁相环波形');

%c6_pllpost;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
kk = 0;
while kk == 0
k = menu('选项',...
        '输入频率和VCO的输出频率',...
        '输入相位和VCO的输出相位',...
        '频率误差图','相位误差图','相位平面图',...
        '相位平面和时域图','退出');
        if k == 1
                plot(t,fin,'r',t,fvco,'k')
                title('输入频率和VCO的输出频率曲线图')
                xlabel('时间/秒');ylabel('频率/Hz');
        elseif k ==2
                pvco=phin-phierror;plot(t,phin,t,pvco)
                title('输入相位和VCO的输出相位曲线图')
                xlabel('时间/秒');ylabel('相位/R');
        elseif k == 3
                plot(t,freqerror);title('频率误差曲线图')
                xlabel('时间/秒');ylabel('频率误差/Hz');
        elseif k == 4
                plot(t,phierror);title('相位误差曲线图')
                xlabel('时间/秒');ylabel('相位差/R');
        elseif k == 5
 %          		 ppplot
            kz = 0;
            while kz == 0
                k = menu('相位平面选项',...
                '延长相位平面',...
                '2pi内的相位平面',...
                '返回');
                if k == 1
                    phierrn = phierror/pi;
                    plot(phierrn,freqerror,'k')
                    title('延长相位平面')
                    xlabel('相位差 /Pi')
                    ylabel('频率差 /Hz')
                    grid

                elseif k == 2
                    pplane(phierror,freqerror,nsettle+1)
                    title('2pi内的相位平面')
                    xlabel('相位差 /Pi')
                    ylabel('频率差 /Hz')
                elseif k == 3
                    kz = 1;
                end
            end
        elseif k == 6   
           		subplot(211);phierrn = phierror/pi;
                plot(phierrn,freqerror,'k');grid;
                title('相位平面曲线图');xlabel('相位差 /Pi');
                ylabel('频率差/Hz');subplot(212)
                plot(t,fin,'k',t,fvco,'k');grid
                title('输入频率和VCO输入频率')
                xlabel('时间/秒');ylabel('频率/Hz');
        elseif k == 7
                kk = 1;
        end
end