1 //
2 // eqlms_cccf_example.c
3 //
4 // Tests least mean-squares (LMS) equalizer (EQ) on a QPSK
5 // signal at the symbol level.
6 //
7
8 #include <stdio.h>
9 #include <stdlib.h>
10 #include <complex.h>
11 #include "liquid.h"
12
13 #define OUTPUT_FILENAME "eqlms_cccf_example.m"
14
main()15 int main() {
16 // options
17 unsigned int n=512; // number of symbols to observe
18 unsigned int ntrain=256;// number of training symbols
19 unsigned int h_len=6; // channel filter length
20 unsigned int p=12; // equalizer order
21
22 // bookkeeping variables
23 float complex d[n]; // data sequence
24 float complex y[n]; // received data sequence (filtered by channel)
25 float complex d_hat[n]; // recovered data sequence
26 float complex h[h_len]; // channel filter coefficients
27 float complex w[p]; // equalizer filter coefficients
28 unsigned int i;
29
30 // create equalizer (default initial coefficients)
31 eqlms_cccf eq = eqlms_cccf_create(NULL,p);
32
33 // create channel filter (random delay taps)
34 h[0] = 1.0f;
35 for (i=1; i<h_len; i++)
36 h[i] = (randnf() + randnf()*_Complex_I) * 0.1f;
37 firfilt_cccf f = firfilt_cccf_create(h,h_len);
38
39 // generate random data signal
40 for (i=0; i<n; i++)
41 d[i] = (rand() % 2 ? 1.0f : -1.0f) +
42 (rand() % 2 ? 1.0f : -1.0f)*_Complex_I;
43
44 // filter data signal through channel
45 for (i=0; i<n; i++) {
46 firfilt_cccf_push(f,d[i]);
47 firfilt_cccf_execute(f,&y[i]);
48 }
49
50 // run equalizer
51 for (i=0; i<p; i++)
52 w[i] = 0;
53 eqlms_cccf_train(eq, w, y, d, ntrain);
54
55 // create filter from equalizer output
56 firfilt_cccf feq = firfilt_cccf_create(w,p);
57
58 // run equalizer filter
59 for (i=0; i<n; i++) {
60 firfilt_cccf_push(feq,y[i]);
61 firfilt_cccf_execute(feq,&d_hat[i]);
62 }
63
64 //
65 // print results
66 //
67 printf("channel:\n");
68 for (i=0; i<h_len; i++)
69 printf(" h(%3u) = %12.8f + j*%12.8f\n", i, crealf(h[i]), cimagf(h[i]));
70
71 printf("equalizer:\n");
72 for (i=0; i<p; i++)
73 printf(" w(%3u) = %12.8f + j*%12.8f\n", i, crealf(w[i]), cimagf(w[i]));
74
75 // compute MSE
76 float complex e;
77 float mse=0.0f;
78 for (i=0; i<n; i++) {
79 // compute mse
80 e = d[i] - d_hat[i];
81 mse += crealf(e*conj(e));
82 }
83 mse /= n;
84 printf("mse: %12.8f\n", mse);
85
86 // clean up objects
87 firfilt_cccf_destroy(f);
88 eqlms_cccf_destroy(eq);
89 firfilt_cccf_destroy(feq);
90
91
92 //
93 // export data to file
94 //
95 FILE * fid = fopen(OUTPUT_FILENAME,"w");
96 fprintf(fid,"%% %s: auto-generated file\n\n", OUTPUT_FILENAME);
97
98 fprintf(fid,"clear all;\n");
99 fprintf(fid,"close all;\n");
100 fprintf(fid,"n=%u;\n",n);
101 fprintf(fid,"ntrain=%u;\n",ntrain);
102 fprintf(fid,"p=%u;\n",p);
103 fprintf(fid,"h_len=%u;\n",h_len);
104
105 for (i=0; i<h_len; i++)
106 fprintf(fid," h(%3u) = %12.4e + j*%12.4e;\n", i+1, crealf(h[i]), cimagf(h[i]));
107
108 for (i=0; i<p; i++)
109 fprintf(fid," w(%3u) = %12.4e + j*%12.4e;\n", i+1, crealf(w[i]), cimagf(w[i]));
110
111 for (i=0; i<n; i++) {
112 fprintf(fid," d(%3u) = %12.4e + j*%12.4e;\n", i+1, crealf(d[i]), cimagf(d[i]));
113 fprintf(fid," y(%3u) = %12.4e + j*%12.4e;\n", i+1, crealf(y[i]), cimagf(y[i]));
114 fprintf(fid," d_hat(%3u) = %12.4e + j*%12.4e;\n", i+1, crealf(d_hat[i]), cimagf(d_hat[i]));
115 }
116
117 // plot results
118 fprintf(fid,"\n\n");
119
120 fprintf(fid,"nfft=512;\n");
121 fprintf(fid,"f=[0:(nfft-1)]/nfft - 0.5;\n");
122 fprintf(fid,"H=20*log10(abs(fftshift(fft(h,nfft))));\n");
123 fprintf(fid,"W=20*log10(abs(fftshift(fft(w,nfft))));\n");
124
125 fprintf(fid,"figure;\n");
126 fprintf(fid,"plot(f,H,'-r',f,W,'-b', f,H+W,'-k','LineWidth',2);\n");
127 fprintf(fid,"xlabel('Normalied Frequency');\n");
128 fprintf(fid,"ylabel('Power Spectral Density [dB]');\n");
129 fprintf(fid,"axis([-0.5 0.5 -10 10]);\n");
130 fprintf(fid,"legend('channel','equalizer','composite',0);\n");
131
132 fprintf(fid,"figure;\n");
133 fprintf(fid,"subplot(2,1,1);\n");
134 fprintf(fid,"hold on;\n");
135 fprintf(fid,"stem(0:(h_len-1),real(h),'-r');\n");
136 fprintf(fid,"stem(0:(p-1), real(w),'-b');\n");
137 fprintf(fid,"hold off;\n");
138 fprintf(fid,"ylabel('Real Coefficients');\n");
139 fprintf(fid,"legend('channel','equalizer',0);\n");
140 fprintf(fid,"axis([-0.25 max(h_len,p)-0.75 -0.5 1.5]);\n");
141 fprintf(fid,"subplot(2,1,2);\n");
142 fprintf(fid,"hold on;\n");
143 fprintf(fid,"stem(0:(h_len-1),imag(h),'-r');\n");
144 fprintf(fid,"stem(0:(p-1), imag(w),'-b');\n");
145 fprintf(fid,"hold off;\n");
146 fprintf(fid,"ylabel('Imag Coefficients');\n");
147 fprintf(fid,"legend('channel','equalizer',0);\n");
148 fprintf(fid,"axis([-0.25 max(h_len,p)-0.75 -0.5 1.5]);\n");
149
150
151 fprintf(fid,"figure;\n");
152 fprintf(fid,"plot(y,'xr',d_hat,'xb');\n");
153 fprintf(fid,"axis('square');\n");
154 fprintf(fid,"xlabel('in-phase');\n");
155 fprintf(fid,"ylabel('quadrature');\n");
156 fprintf(fid,"legend('received','equalized',1');\n");
157
158 fclose(fid);
159 printf("results written to %s.\n",OUTPUT_FILENAME);
160
161 return 0;
162 }
163