1 //
2 // iirdes_pll_example.c
3 //
4 // This example demonstrates 2nd-order IIR phase-locked loop filter
5 // design with a practical simulation.
6 //
7 // SEE ALSO: nco_pll_example.c
8 // nco_pll_modem_example.c
9 //
10
11 #include <stdio.h>
12 #include <stdlib.h>
13 #include <time.h>
14 #include <math.h>
15 #include <getopt.h>
16
17 #include "liquid.h"
18
19 #define OUTPUT_FILENAME "iirdes_pll_example.m"
20
21 // print usage/help message
usage()22 void usage()
23 {
24 printf("iirdes_pll_example [options]\n");
25 printf(" u/h : print usage\n");
26 printf(" b : pll bandwidth, default: 0.01\n");
27 printf(" z : zeta (damping factor), default 0.70711\n");
28 printf(" K : loop filter gain, default: 1000\n");
29 printf(" n : number of samples, default: 512\n");
30 printf(" p : phase offset (radians), default: pi/4\n");
31 printf(" f : frequency offset (radians), default: 0.3\n");
32 }
33
main(int argc,char * argv[])34 int main(int argc, char*argv[]) {
35 srand( time(NULL) );
36
37 // options
38 float phase_offset = M_PI / 4.0f; // phase offset
39 float frequency_offset = 0.3f; // frequency offset
40 float pll_bandwidth = 0.01f; // PLL bandwidth
41 float zeta = 1/sqrtf(2.0f); // PLL damping factor
42 float K = 1000.0f; // PLL loop gain
43 unsigned int n=512; // number of iterations
44
45 int dopt;
46 while ((dopt = getopt(argc,argv,"uhb:z:K:n:p:f:")) != EOF) {
47 switch (dopt) {
48 case 'u':
49 case 'h': usage(); return 0;
50 case 'b': pll_bandwidth = atof(optarg); break;
51 case 'z': zeta = atof(optarg); break;
52 case 'K': K = atof(optarg); break;
53 case 'n': n = atoi(optarg); break;
54 case 'p': phase_offset = atof(optarg); break;
55 case 'f': frequency_offset= atof(optarg); break;
56 default:
57 exit(1);
58 }
59 }
60 unsigned int d=n/32; // print every "d" lines
61
62 // validate input
63 if (pll_bandwidth <= 0.0f) {
64 fprintf(stderr,"error: bandwidth must be greater than 0\n");
65 exit(1);
66 } else if (zeta <= 0.0f) {
67 fprintf(stderr,"error: damping factor must be greater than 0\n");
68 exit(1);
69 } else if (K <= 0.0f) {
70 fprintf(stderr,"error: loop gain must be greater than 0\n");
71 exit(1);
72 }
73
74 // data arrays
75 float complex x[n]; // input complex sinusoid
76 float complex y[n]; // output complex sinusoid
77 float phase_error[n]; // output phase error
78
79 // generate PLL filter
80 float b[3];
81 float a[3];
82 iirdes_pll_active_lag(pll_bandwidth, zeta, K, b, a);
83 iirfilt_rrrf pll = iirfilt_rrrf_create(b,3,a,3);
84 iirfilt_rrrf_print(pll);
85
86 unsigned int i;
87 float phi;
88 for (i=0; i<n; i++) {
89 phi = phase_offset + i*frequency_offset;
90 x[i] = cexpf(_Complex_I*phi);
91 }
92
93 // run loop
94 float theta = 0.0f;
95 y[0] = 1.0f;
96 for (i=0; i<n; i++) {
97
98 // generate complex sinusoid
99 y[i] = cexpf(_Complex_I*theta);
100
101 // compute phase error
102 phase_error[i] = cargf(x[i]*conjf(y[i]));
103
104 // update pll
105 iirfilt_rrrf_execute(pll, phase_error[i], &theta);
106
107 // print phase error
108 if ((i)%d == 0 || i==n-1 || i==0)
109 printf("%4u : phase error = %12.8f\n", i, phase_error[i]);
110 }
111
112 // destroy filter object
113 iirfilt_rrrf_destroy(pll);
114
115 // write output file
116 FILE * fid = fopen(OUTPUT_FILENAME,"w");
117 fprintf(fid,"%% %s : auto-generated file\n", OUTPUT_FILENAME);
118 fprintf(fid,"clear all;\n");
119 fprintf(fid,"close all;\n");
120 fprintf(fid,"n = %u;\n", n);
121 fprintf(fid,"x = zeros(1,n);\n");
122 fprintf(fid,"y = zeros(1,n);\n");
123 for (i=0; i<n; i++) {
124 fprintf(fid,"x(%4u) = %12.4e + j*%12.4e;\n", i+1, crealf(x[i]), cimagf(x[i]));
125 fprintf(fid,"y(%4u) = %12.4e + j*%12.4e;\n", i+1, crealf(y[i]), cimagf(y[i]));
126 fprintf(fid,"e(%4u) = %12.4e;\n", i+1, phase_error[i]);
127 }
128 fprintf(fid,"t=0:(n-1);\n");
129 fprintf(fid,"figure;\n");
130 fprintf(fid,"subplot(2,1,1);\n");
131 fprintf(fid," plot(t,real(x),t,real(y));\n");
132 fprintf(fid," xlabel('time');\n");
133 fprintf(fid," ylabel('real');\n");
134 fprintf(fid,"subplot(2,1,2);\n");
135 fprintf(fid," plot(t,imag(x),t,imag(y));\n");
136 fprintf(fid," xlabel('time');\n");
137 fprintf(fid," ylabel('imag');\n");
138
139 fprintf(fid,"figure;\n");
140 fprintf(fid,"plot(t,e);\n");
141 fprintf(fid,"xlabel('time');\n");
142 fprintf(fid,"ylabel('phase error');\n");
143 fprintf(fid,"grid on;\n");
144
145 fclose(fid);
146 printf("results written to %s.\n",OUTPUT_FILENAME);
147
148 printf("done.\n");
149 return 0;
150 }
151