1 //
2 // eqlms_cccf_test.c
3 //
4 // Tests least mean-squares (LMS) equalizer (EQ).
5 //
6
7 #include <stdio.h>
8 #include <stdlib.h>
9 #include <string.h>
10 #include <math.h>
11 #include <complex.h>
12 #include "liquid.h"
13
14 #define OUTPUT_FILENAME "eqlms_cccf_test.m"
15
main()16 int main() {
17 // options
18 unsigned int k=2; // samples/symbol
19 unsigned int m=3;
20 float beta = 0.3f;
21 //float dt = 0.4f; // timing offset
22 unsigned int num_symbols = 256;
23 unsigned int hc_len=1; // channel impulse response length
24 unsigned int p=13; // equalizer filter length
25 float mu = 0.05f;
26
27 // TODO : validate input
28
29 // derived values
30 unsigned int num_samples = k*num_symbols;
31
32 // data arrays
33 float complex s[num_symbols]; // original QPSK symbols
34 float complex x[num_samples]; // interpolated samples
35 float complex y[num_samples]; // received samples
36 float complex z[num_samples]; // recovered samples
37
38 // generate data sequence
39 unsigned int i;
40 for (i=0; i<num_symbols; i++) {
41 s[i] = (rand() % 2 ? 1.0f : -1.0f) +
42 (rand() % 2 ? 1.0f : -1.0f) * _Complex_I;
43 }
44
45 // interpolate
46 unsigned int h_len = 2*k*m+1;
47 float hm[h_len];
48 #if 1
49 // design GMSK
50 float c0 = 1.0f / sqrtf(logf(2.0f));
51 for (i=0; i<h_len; i++) {
52 float t = (float)(i) / (float)(k) - (float)(m);
53 //float sig = 0.7f;
54 //hm[i] = expf(-t*t / (2*sig*sig) );
55
56 hm[i] = liquid_Qf(2*M_PI*beta*(t-0.5f)*c0) -
57 liquid_Qf(2*M_PI*beta*(t+0.5f)*c0);
58 printf("h(%3u) = %12.8f\n", i+1, hm[i]);
59 }
60 #else
61 design_rrc_filter(k,m,beta,dt,hm);
62 #endif
63 firinterp_crcf interp = firinterp_crcf_create(k,hm,h_len);
64 for (i=0; i<num_symbols; i++)
65 firinterp_crcf_execute(interp, s[i], &x[i*k]);
66 firinterp_crcf_destroy(interp);
67
68 // generate channel filter
69 float complex hc[hc_len];
70 for (i=0; i<hc_len; i++)
71 hc[i] = (i==0) ? 1.0f : 0.1f*randnf()*cexpf(_Complex_I*2*M_PI*randf());
72
73 // push signal through channel...
74 firfilt_cccf fchannel = firfilt_cccf_create(hc,hc_len);
75 for (i=0; i<num_samples; i++) {
76 firfilt_cccf_push(fchannel, x[i]);
77 firfilt_cccf_execute(fchannel, &y[i]);
78 }
79 firfilt_cccf_destroy(fchannel);
80
81 // initialize equalizer
82 #if 1
83 // initialize with delta at center
84 float complex h[p]; // coefficients
85 unsigned int p0 = (p-1)/2; // center index
86 for (i=0; i<p; i++)
87 h[i] = (i==p0) ? 1.0f : 0.0f;
88 #else
89 // initialize with square-root raised cosine
90 p = 2*k*m+1;
91 float h_tmp[p];
92 float complex h[p];
93 design_rrc_filter(k,m,beta,0.0f,h_tmp);
94 for (i=0; i<p; i++)
95 h[i] = h_tmp[i] / k;
96 #endif
97
98 // run equalizer
99 float complex w0[p];
100 float complex w1[p];
101 memmove(w0, h, p*sizeof(float complex));
102 windowcf buffer = windowcf_create(p);
103 for (i=0; i<num_samples; i++) {
104 // push value into buffer
105 windowcf_push(buffer, y[i]);
106
107 // compute d_hat
108 float complex d_hat = 0.0f;
109 float complex * r;
110 windowcf_read(buffer, &r);
111 unsigned int j;
112 for (j=0; j<p; j++)
113 d_hat += conjf(w0[j]) * r[j];
114
115 // store in output
116 z[i] = d_hat;
117
118 // check to see if buffer is full, return if not
119 if (i <= p) continue;
120
121 // decimate by k
122 if ( (i%k) != 0 ) continue;
123
124 // estimate transmitted QPSK symbol
125 float complex d_prime = (crealf(d_hat) > 0.0f ? 1.0f : -1.0f) +
126 (cimagf(d_hat) > 0.0f ? 1.0f : -1.0f) * _Complex_I;
127
128 // compute error
129 float complex alpha = d_prime - d_hat;
130
131 // compute signal energy
132 float ex2 = 1.414f;
133
134 // update weighting vector
135 for (j=0; j<p; j++)
136 w1[j] = w0[j] + mu*conjf(alpha)*r[j]/ex2;
137
138 // copy new filter values
139 memmove(w0, w1, p*sizeof(complex float));
140 }
141
142 // destroy additional objects
143 windowcf_destroy(buffer);
144
145 // print results to file
146 FILE * fid = fopen(OUTPUT_FILENAME,"w");
147 fprintf(fid,"%% %s : auto-generated file\n\n", OUTPUT_FILENAME);
148 fprintf(fid,"clear all\n");
149 fprintf(fid,"close all\n");
150
151 fprintf(fid,"k = %u;\n", k);
152 fprintf(fid,"m = %u;\n", m);
153 fprintf(fid,"num_symbols = %u;\n", num_symbols);
154 fprintf(fid,"num_samples = num_symbols*k;\n");
155
156 // save samples
157 fprintf(fid,"x = zeros(1,num_samples);\n");
158 fprintf(fid,"y = zeros(1,num_samples);\n");
159 fprintf(fid,"z = zeros(1,num_samples);\n");
160 for (i=0; i<num_samples; i++) {
161 fprintf(fid,"x(%3u) = %12.4e + j*%12.4e;\n", i+1, crealf(x[i]), cimagf(x[i]));
162 fprintf(fid,"y(%3u) = %12.4e + j*%12.4e;\n", i+1, crealf(y[i]), cimagf(y[i]));
163 fprintf(fid,"z(%3u) = %12.4e + j*%12.4e;\n", i+1, crealf(z[i]), cimagf(z[i]));
164 }
165 fprintf(fid,"t = 1:num_samples;\n");
166 fprintf(fid,"tsym = 1:k:num_samples;\n");
167 fprintf(fid,"figure;\n");
168 fprintf(fid,"subplot(2,1,1), plot(t, real(z), tsym, real(z(tsym)),'x');\n");
169 fprintf(fid,"subplot(2,1,2), plot(t, imag(z), tsym, imag(z(tsym)),'x');\n");
170
171 // plot symbols
172 fprintf(fid,"tsym0 = tsym([1:length(tsym)/2]);\n");
173 fprintf(fid,"tsym1 = tsym((length(tsym)/2):end);\n");
174 fprintf(fid,"figure;\n");
175 fprintf(fid,"plot(real(z(tsym0)),imag(z(tsym0)),'x','Color',[0.7 0.7 0.7],...\n");
176 fprintf(fid," real(z(tsym1)),imag(z(tsym1)),'x','Color',[0.0 0.0 0.0]);\n");
177 fprintf(fid,"xlabel('In-phase');\n");
178 fprintf(fid,"ylabel('Quadrature');\n");
179 fprintf(fid,"axis([-1 1 -1 1]*1.5);\n");
180 fprintf(fid,"axis square;\n");
181 fprintf(fid,"grid on;\n");
182
183 // save filter coefficients
184 for (i=0; i<h_len; i++) fprintf(fid,"hm(%3u) = %12.4e;\n", i+1, hm[i]);
185 for (i=0; i<p; i++) fprintf(fid,"hp(%3u) = %12.4e + j*%12.4e;\n", i+1, crealf(w0[i]), cimagf(w0[i]));
186 for (i=0; i<hc_len; i++)fprintf(fid,"hc(%3u) = %12.4e + j*%12.4e;\n", i+1, crealf(hc[i]), cimagf(hc[i]));
187
188 fprintf(fid,"nfft = 1024;\n");
189 fprintf(fid,"f = [0:(nfft-1)]/nfft - 0.5;\n");
190 fprintf(fid,"Hm = 20*log10(abs(fftshift(fft(hm/k,nfft))));\n");
191 fprintf(fid,"Hc = 20*log10(abs(fftshift(fft(hc,nfft))));\n");
192 fprintf(fid,"Hp = 20*log10(abs(fftshift(fft(hp,nfft))));\n");
193 fprintf(fid,"G = Hm + Hp + Hc;\n");
194 fprintf(fid,"figure;\n");
195 fprintf(fid,"plot(f,Hm, f,Hc, f,Hp, f,G,'-k','LineWidth',2);\n");
196 fprintf(fid,"xlabel('Normalized Frequency');\n");
197 fprintf(fid,"ylabel('Power Spectral Density');\n");
198 fprintf(fid,"legend('transmit','channel','equalizer','total');\n");
199 fprintf(fid,"axis([-0.5 0.5 -10 10]);\n");
200 fprintf(fid,"grid on;\n");
201
202 fclose(fid);
203
204 printf("results written to '%s'\n", OUTPUT_FILENAME);
205
206 printf("done.\n");
207 return 0;
208 }
209