%%DD算法---频率捕获范围
%仿真参数
N = 500000;
M = 256;
SNR = 30;
samplingFreq = 5000000;
carrFreqOffset = -300000:20000:300000;
carrPhsOffset = 0;
h1 = modem.qammod('M',2^8, 'SymbolOrder', 'Gray');
h2 = modem.qamdemod('M', 2^8, 'SymbolOrder', 'Gray');

%鉴相器参量
DifferOfPha = zeros(1,N);
DD_bitsOutput = zeros(1,N);
DD_DifferOfPha = zeros(1,N);
Z = zeros(1,N);
Y = zeros(1,256);

%锁定检测器参量
lamuda = 0.7;
beita = 0.6;
Ncounter = 0;
Track_sign = 0;
MeanOfY = 0;
Lock_N = zeros(1,length(carrFreqOffset));

%环路滤波器及NCO参量
fs = samplingFreq;
fn = 50000;
yita = 0.5;
wn = 2*pi*fn/fs;
Kp = 2*yita*wn;
Ki = wn^2;
PraZ = 0;
PhasControl = 0;
PhaseOfNCO = 0;

%环路捕获频率
PreAcqFreq = 0;
RealAcqFreq = zeros(1,length(carrFreqOffset));

%通信过程仿真
for fre = 1:1:length(carrFreqOffset)
bitSrc = randi([0 M-1],1,N);
bitsTransmit = modulate(h1,bitSrc);
phaseStep = carrFreqOffset(fre) / samplingFreq;
phaseVar = phaseStep * (0:1:length(bitsTransmit)-1);
aftFreOffset = bitsTransmit .* exp(1j*(2*pi*phaseVar+carrPhsOffset));
bitsnoise = awgn(aftFreOffset,SNR,'measured');
for m=1:N
%%PD
DifferOfPha(m) = bitsnoise(m)*exp(-1j*PhaseOfNCO);
DD_bitsOutput(m) = demodulate(h2,DifferOfPha(m));
DD_DifferOfPha(m) = modulate(h1,DD_bitsOutput(m));

Ncounter=Ncounter+1;
if(abs(DD_DifferOfPha(m)-DifferOfPha(m))Y(Ncounter)=1;
else
Y(Ncounter)=0;
end
if(Ncounter==256)
MeanOfY = mean(Y);
Ncounter = 0;
end
if(Track_sign==0)
Z(m) = imag(DifferOfPha(m)/DD_DifferOfPha(m));
if(MeanOfY>yita)
Track_sign = 1;
Lock_N(fre) = m;
end
else
Z(m) = imag(DifferOfPha(m)/DD_DifferOfPha(m));
end

%%Loop Filter
Phaz = Kp*Z(m) + PhasControl;
PhasControl = Ki*Z(m) + PhasControl;

%%NCO
PhaseOfNCO = PhaseOfNCO + Phaz;

%%acqucisition frequency
Acqfreq = 0.01*PhasControl + 0.99*PreAcqFreq;
PreAcqFreq = Acqfreq;
RealAcqFreq(fre) = PreAcqFreq/2/pi*5000000;
end
end

%figure
figure(1);
i1 = 1:1:length(carrFreqOffset);
plot((i1-16)*20000,RealAcqFreq,'k-*','linewidth',2);
xlabel('实际频率偏移/Hz');
ylabel('环路捕获频率/Hz');
grid on;

if 0
figure(2);
i2 = 1:1:length(carrFreqOffset);
plot((i1-16)*20000,Lock_N,'k-*','linewidth',2);
xlabel('实际频率偏移/Hz');
ylabel('环路工作时间/T');
end