From: HardySpicer on
On Apr 16, 10:31 am, "smc123" <xmacleod(a)n_o_s_p_a_m.yahoo.com> wrote:
> I'm looking for guidance on writing a software PLL for a signal acquired
> from a data acquisition board:
> The daq board will be sampling at 10kSPS.
> The signals to which I'd like to lock on to are from 5Hz to 1kHz, sine or
> square wave.
> The output of the PLL will be used to multiply a signal acquired on an
> adjacent channel of the daq board.  Lock-in amplifier application.
> The processing of the signals will be done on a PC-based software app.
>
> Thanks for your help

%Phase Locked Loop

%This simulates an analogue PLL
%The VCO is sinusoidal but can be made a square wave
%Suggested parameters: Sample at 10000Hz,Carrier freq 1000Hz = VCO
free running freq
% Make baseband freq 1Hz and frequency deviation 10Hz.
% Make run time 1 second
%When you run this you will see the demodulated sine wave (1Hz) as
figure 1.
% It will have 2fcarrier (2fc) superimposed on top of it (as in a real
PLL)

%First Generate Test Signal
%Frequency Modulation of a sine wave is simplest!
fs=input('Input Sampling Frequency (Hz)');
fb=input('Baseband Frequency (FM)');
fc=input('Carrier Frequency (FM)');
fset=input('Free running frequency for VCO(make this normally =
carrier freq)');
delf=input('Frequency Deviation (Hz)');
%Compute FM Modulation Index
beta=delf/fb;
tmax=input('Run Time (Secs)');
npoints=tmax*fs;

%%%%%%%%%%%%%%%%%%%%%%%%%%
%Some arrays inialised
pout=zeros(npoints,1);
fm=zeros(npoints,1);
vco=zeros(npoints,1);
ff1=zeros(npoints,1);
ff2=zeros(npoints,1);
%%%%%%%%%%%%%%%%%%%%%%%%%

%Generate FM

for i=1:1:npoints
fm(i)=cos(2*pi*i*(fc/fs) + beta* sin(2*pi*i*(fb/fs) ) );
end


%Unity Gain Crossover Frequency (Hz) for Bode Plot
fcp=(2*fc)/10;

%Make span ratio 10:1 gives a phase margin of around 57 degrees
%Change this and you change stability phase margin is arcsin[(span-1)/
(span+1)]
spanp=10;
%This function computes phase advance parameters
[ap,bp,delp]=plead(spanp,fs,fcp);

%Compute Overall Gain

gainp=2*(1-cos (2*pi*(fcp/fs) ) )/sqrt(spanp);
%Divide up between Integrators
gainip=sqrt(gainp);

%Initialise just in case..
p0=0.0;
p1=0.0;
q0=0.0;
q1=0.0;
r0=0.0;
r1=0.0;
fu0=0.0;


%Main Loop for Phase Locked Loop

for i=1:1:npoints

%Controller (Filter)
%First Integrator
p1=p0;
p0=p1+fu0*gainip;

%fu0 is Phase Detector Output

%Second Integrator
q1=q0;
q0=q1+gainip*p0;

%Phase Advance
r1=r0;
r0=-r1*bp + delp*(q0+ap*q1);


%VCO Output
vco(i)=cos( (2*pi*i*(fset/fs) ) + r0 ) ;

%Now Limit VCO Output to +1.0 or - 1.0
%For square wave VCO uncomment the following
%if(vco(i)>0)vco(i)=1.0;end
%if(vco(i)<0)vco(i)=-1.0;end









%Phase Detector (A simple multiplier)
fu0=vco(i)*fm(i);


%PLL Output is first Integrator Output
pout(i)=p0;


end

figure(1)

plot(pout)


function [a,b,del]=plead(span,fs,fc)
%Phase Lead
% Returns zero coefficient a
% Pole coefficient b
% and gain del
%span ratio is span
%fs is sampling frequency
%fc is frenquency for maximum phase advance
thetac=2*pi*(fc/fs);
sqs=sqrt(span);


x=cos(thetac);

y=sin(thetac);

num1=sqs*(1+span)*y-2*span;

den1=sqs*(1-span)*y + 2*span*x;


a=num1/den1;

beta=(1-a)/(1+a);

b=(span-beta)/(span+beta);

del=(1+b)/(1+a);

%Th Thats all folks....