锁频环FLL的matlab仿真

%---------------------------------------------------------------------%
%-------------------------FLL-----------------------------------------%
%---------------------------------------------------------------------%

clc;

    fs = 16.368e6;				%[Hz]	采样频率
	fi = 4.090e6;				%[Hz]	中频频率
    T = 0.001;					%[s]	积分清除时间1ms
    n_accum = T * fs;			%     	积分清除点数
    CNR = 44;					%	  	载噪比44dB==-130dBmw
    SNR = CNR -10*log10（fs）;	%	  	由载噪比得到信噪比
    addnoise = 1;				% 1: 加噪声; 0: 不加噪声
    fd0 = 300;					%[Hz] 	多普勒初始频率
    detfd = 1000;				%[Hz/s]	多普勒变化率
    kf = 2^32/fs;				%		频率控制字转换系数
    time = 0.5;					%		仿真时间和仿真次数
    loopnum = fix（time/T）;		%		仿真点数/环路更新次数
    b_fll = 0.9;				%[Hz]	锁频环带宽
    wn = 1.89 * b_fll;			%		阻尼系数为0.707下中心频率
    gain_nco = 2*pi*fs/2^32;	%		NCO增益
    gain_disctim = 1;			%		鉴相器增益
    K = gain_nco * gain_disctim;%		环路增益：NCO*鉴相器    

    k1 = 0.3;	% 修改自己的滤波系数
    k2 = 0.2;				%修改自己的滤波系数
	
	% New FLL-------------------------------------------------------------
	fllDampingRatio		= 0.7;		% 
	fllNoiseBandwidth	= 3;		%[Hz]
	fllCoherentT		= 0.001;	% Coherent Time
	% FLL parameter calculation accroding to settings-------
	[tau1Fll tau2Fll] = calcLoopCoef（fllNoiseBandwidth ...
									fllDampingRatio K）;
	oldFllNco   = 0.0;
	oldFllError = 0.0;
	oldoldFllNco = 0.0;

    endlasttime = 0;
    
    fd_estimate = 0;
    phaseadder = 0;
    
    past_I_accum = 0;
    past_Q_accum = 0;
    
    freqoutdlf = 0;
    past_past_freqoutdlf = 0;
    past_freqoutdlf = 0;
    past_freqerror = 0;
	
	freq = zeros（1loopnum）;

for index = 1 : loopnum

	    t = endlasttime + （1:n_accum）/fs;
            phase = 2*pi* [（fi+fd0）*t + 1/2*detfd*t.^2];
%            phase = 2*pi*（fi+fd）*（endlasttime + （1:n_accum）/fs）;
            endlasttime = endlasttime + T;
            freq（index） = fd0 + detfd*endlasttime;
            cos_signal = cos（phase）;
        if addnoise == 0
            cos_source = cos_signal;
        e