costa.m

锁相环的仿真以及含有文本

clc;
close all;
clear all;

file='BPSK_data.dat';
fid=fopen(file,'r');
sig=fread(fid,'float32');
fclose(fid); 



%%

U10 = 1; % Input Reference signal amplitude
U20 = 1; % Output NCO signal amplitude
F0 = 61; % Center Frequency of NCO
Fs = 200 % Sampling freqency
Ts = 1/Fs; % Sampling duration


%%  画出原始信号及其频谱，观测形态是否正常 

x=sig(1:2048);
t=0:Fs/length(x):Fs-Fs/length(x);
figure;
plot(x);
title('待解调信号时域波形');
sp=20*log10(abs(fft(x)));
figure;
plot(t,sp);
title('待解调信号频谱');


%%
lpCoeff=importdata('loCoeff.txt');
lpCoeff=lpCoeff';
firOrder=length(lpCoeff);
waitFilBufI=zeros(firOrder,1);  
waitFilBufQ=zeros(firOrder,1); 


% Calculation of PLL parameters
zeta = 0.707; % Damping factor
Wl = 30*pi; % Lock range defined according to expected
% input characteristics
Wn = Wl/2/zeta % Natural frequency calculated from lock range
Kd = 1; % Phase detector gain
Ko = 1; % Loop filter gain

% Loop filter constants and lead-lag gain constants calculated
tau1 = Ko*Kd/(Wn^2);
tau2 = 2*zeta/Wn;
a1 = -1;
b0 = (Ts/(2*tau1))*(1+1/(tan(Ts/(2*tau2)))); % Lead constant gain
b1 = (Ts/(2*tau1))*(1-1/(tan(Ts/(2*tau2)))); % Lag constant gain


% Initialization for first cycle of feedback loop
phi = 0; % Output phase of the NCO
theta = 0; % Phase error of the two signals
Uf=0; % Output signal of the loop filter
Imix=0;
Qmix=0;
mix=0;
theta=0;
F1 = 2*pi*F0 % Output frequency of NCO

% General case phase lock loop

for n = 2:1:length(sig)
  
  phi(n) = phi(n-1)+F1(n-1)*Ts; % Phase error from current cycle
  
  loI=cos(phi(n)); %NCO同相输出
  loQ=sin(phi(n)); %NCO正交输出

  Imix(n)=sig(n)*loI;
  Qmix(n)=sig(n)*loQ;    
  
  %fir低通滤波  
  for j=firOrder:-1:2
      waitFilBufI(j)=waitFilBufI(j-1);
  end;
  waitFilBufI(1)=Imix(n);    
  Ilp(n)=lpCoeff*waitFilBufI;

  for j=firOrder:-1:2
      waitFilBufQ(j)=waitFilBufQ(j-1);
  end;
  waitFilBufQ(1)=Qmix(n);
  Qlp(n)=lpCoeff*waitFilBufQ;

  theta(n)=Ilp(n)*Qlp(n);    %相位误差
  
  Uf(n) = -a1*Uf(n-1)+b0*theta(n)+b1*theta(n-1);        
  
  F1(n) = 2*pi*F0+Uf(n);
       
end 

figure;
plot(theta);

figure;
plot(Uf);