QAM盲载波恢复算法---DD算法(频率捕获范围代码)

%%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