1 //
2 // eqlms_cccf_blind_example.c
3 //
4 // This example tests the least mean-squares (LMS) equalizer (EQ) on a
5 // signal with an unknown modulation and carrier frequency offset. That
6 // is, the equalization is done completely blind of the modulation
7 // scheme or its underlying data set. The error estimate assumes a
8 // constant modulus linear modulation scheme. This works surprisingly
9 // well even more amplitude-modulated signals, e.g. 'qam16'.
10 //
11
12 #include <stdio.h>
13 #include <stdlib.h>
14 #include <string.h>
15 #include <math.h>
16 #include <complex.h>
17 #include <getopt.h>
18 #include <time.h>
19 #include "liquid.h"
20
21 #define OUTPUT_FILENAME "eqlms_cccf_blind_example.m"
22
23 // print usage/help message
usage()24 void usage()
25 {
26 printf("Usage: eqlms_cccf_blind_example [OPTION]\n");
27 printf(" h : print help\n");
28 printf(" n : number of symbols, default: 500\n");
29 printf(" s : SNR [dB], default: 30\n");
30 printf(" c : number of channel filter taps (minimum: 1), default: 5\n");
31 printf(" k : samples/symbol, default: 2\n");
32 printf(" m : filter semi-length (symbols), default: 4\n");
33 printf(" b : filter excess bandwidth factor, default: 0.3\n");
34 printf(" p : equalizer semi-length (symbols), default: 3\n");
35 printf(" u : equalizer learning rate, default; 0.05\n");
36 printf(" M : modulation scheme (qpsk default)\n");
37 liquid_print_modulation_schemes();
38 }
39
main(int argc,char * argv[])40 int main(int argc, char*argv[])
41 {
42 //srand(time(NULL));
43
44 // options
45 unsigned int num_symbols=800; // number of symbols to observe
46 float SNRdB = 30.0f; // signal-to-noise ratio [dB]
47 float fc = 0.002f; // carrier offset
48 unsigned int hc_len=5; // channel filter length
49 unsigned int k=2; // matched filter samples/symbol
50 unsigned int m=3; // matched filter delay (symbols)
51 float beta=0.3f; // matched filter excess bandwidth factor
52 unsigned int p=3; // equalizer length (symbols, hp_len = 2*k*p+1)
53 float mu = 0.08f; // equalizer learning rate
54
55 // modulation type/depth
56 modulation_scheme ms = LIQUID_MODEM_QPSK;
57
58 int dopt;
59 while ((dopt = getopt(argc,argv,"hn:s:c:k:m:b:p:u:M:")) != EOF) {
60 switch (dopt) {
61 case 'h': usage(); return 0;
62 case 'n': num_symbols = atoi(optarg); break;
63 case 's': SNRdB = atof(optarg); break;
64 case 'c': hc_len = atoi(optarg); break;
65 case 'k': k = atoi(optarg); break;
66 case 'm': m = atoi(optarg); break;
67 case 'b': beta = atof(optarg); break;
68 case 'p': p = atoi(optarg); break;
69 case 'u': mu = atof(optarg); break;
70 case 'M':
71 ms = liquid_getopt_str2mod(optarg);
72 if (ms == LIQUID_MODEM_UNKNOWN) {
73 fprintf(stderr,"error: %s, unknown/unsupported modulation scheme '%s'\n", argv[0], optarg);
74 return 1;
75 }
76 break;
77 default:
78 exit(1);
79 }
80 }
81
82 // validate input
83 if (num_symbols == 0) {
84 fprintf(stderr,"error: %s, number of symbols must be greater than zero\n", argv[0]);
85 exit(1);
86 } else if (hc_len == 0) {
87 fprintf(stderr,"error: %s, channel must have at least 1 tap\n", argv[0]);
88 exit(1);
89 } else if (k < 2) {
90 fprintf(stderr,"error: %s, samples/symbol must be at least 2\n", argv[0]);
91 exit(1);
92 } else if (m == 0) {
93 fprintf(stderr,"error: %s, filter semi-length must be at least 1 symbol\n", argv[0]);
94 exit(1);
95 } else if (beta < 0.0f || beta > 1.0f) {
96 fprintf(stderr,"error: %s, filter excess bandwidth must be in [0,1]\n", argv[0]);
97 exit(1);
98 } else if (p == 0) {
99 fprintf(stderr,"error: %s, equalizer semi-length must be at least 1 symbol\n", argv[0]);
100 exit(1);
101 } else if (mu < 0.0f || mu > 1.0f) {
102 fprintf(stderr,"error: %s, equalizer learning rate must be in [0,1]\n", argv[0]);
103 exit(1);
104 }
105
106 // derived values
107 unsigned int hm_len = 2*k*m+1; // matched filter length
108 unsigned int hp_len = 2*k*p+1; // equalizer filter length
109 unsigned int num_samples = k*num_symbols;
110
111 // bookkeeping variables
112 float complex syms_tx[num_symbols]; // transmitted data symbols
113 float complex x[num_samples]; // interpolated time series
114 float complex y[num_samples]; // channel output
115 float complex z[num_samples]; // equalized output
116 float complex syms_rx[num_symbols]; // received data symbols
117
118 float hm[hm_len]; // matched filter response
119 float complex hc[hc_len]; // channel filter coefficients
120 float complex hp[hp_len]; // equalizer filter coefficients
121
122 unsigned int i;
123
124 // generate matched filter response
125 liquid_firdes_prototype(LIQUID_FIRFILT_RRC, k, m, beta, 0.0f, hm);
126 firinterp_crcf interp = firinterp_crcf_create(k, hm, hm_len);
127
128 // create the modem objects
129 modem mod = modem_create(ms);
130 modem demod = modem_create(ms);
131 unsigned int M = 1 << modem_get_bps(mod);
132
133 // generate channel impulse response, filter
134 hc[0] = 1.0f;
135 for (i=1; i<hc_len; i++)
136 hc[i] = 0.09f*(randnf() + randnf()*_Complex_I);
137 firfilt_cccf fchannel = firfilt_cccf_create(hc, hc_len);
138
139 // generate random symbols
140 for (i=0; i<num_symbols; i++)
141 modem_modulate(mod, rand()%M, &syms_tx[i]);
142
143 // interpolate
144 for (i=0; i<num_symbols; i++)
145 firinterp_crcf_execute(interp, syms_tx[i], &x[i*k]);
146
147 // push through channel
148 float nstd = powf(10.0f, -SNRdB/20.0f);
149 for (i=0; i<num_samples; i++) {
150 firfilt_cccf_push(fchannel, x[i]);
151 firfilt_cccf_execute(fchannel, &y[i]);
152
153 // add carrier offset and noise
154 y[i] *= cexpf(_Complex_I*2*M_PI*fc*i);
155 y[i] += nstd*(randnf() + randnf()*_Complex_I)*M_SQRT1_2;
156 }
157
158 // push through equalizer
159 // create equalizer, intialized with square-root Nyquist filter
160 eqlms_cccf eq = eqlms_cccf_create_rnyquist(LIQUID_FIRFILT_RRC, k, p, beta, 0.0f);
161 eqlms_cccf_set_bw(eq, mu);
162
163 // get initialized weights
164 eqlms_cccf_get_weights(eq, hp);
165
166 // filtered error vector magnitude (emperical RMS error)
167 float evm_hat = 0.03f;
168
169 // nco/pll for phase recovery
170 nco_crcf nco = nco_crcf_create(LIQUID_VCO);
171 nco_crcf_pll_set_bandwidth(nco, 0.02f);
172
173 float complex d_hat = 0.0f;
174 unsigned int num_symbols_rx = 0;
175 for (i=0; i<num_samples; i++) {
176 // print filtered evm (emperical rms error)
177 if ( ((i+1)%50)==0 )
178 printf("%4u : rms error = %12.8f dB\n", i+1, 10*log10(evm_hat));
179
180 eqlms_cccf_push(eq, y[i]);
181 eqlms_cccf_execute(eq, &d_hat);
182
183 // store output
184 z[i] = d_hat;
185
186 // decimate by k
187 if ( (i%k) != 0 ) continue;
188
189 // update equalizer independent of the signal: estimate error
190 // assuming constant modulus signal
191 eqlms_cccf_step_blind(eq, d_hat);
192
193 // apply carrier recovery
194 float complex v;
195 nco_crcf_mix_down(nco, d_hat, &v);
196
197 // save resulting data symbol
198 if (num_symbols_rx < num_symbols)
199 syms_rx[num_symbols_rx++] = v;
200
201 // demodulate
202 unsigned int sym_out; // output symbol
203 float complex d_prime; // estimated input sample
204 modem_demodulate(demod, v, &sym_out);
205 modem_get_demodulator_sample(demod, &d_prime);
206 float phase_error = modem_get_demodulator_phase_error(demod);
207
208 // update pll
209 nco_crcf_pll_step(nco, phase_error);
210
211 // update rx nco object
212 nco_crcf_step(nco);
213
214 // update filtered evm estimate
215 float evm = crealf( (d_prime-v)*conjf(d_prime-v) );
216 evm_hat = 0.98f*evm_hat + 0.02f*evm;
217 }
218
219 // get equalizer weights
220 eqlms_cccf_get_weights(eq, hp);
221
222 // destroy objects
223 eqlms_cccf_destroy(eq);
224 nco_crcf_destroy(nco);
225 firinterp_crcf_destroy(interp);
226 firfilt_cccf_destroy(fchannel);
227 modem_destroy(mod);
228 modem_destroy(demod);
229
230 //
231 // export output
232 //
233 FILE * fid = fopen(OUTPUT_FILENAME,"w");
234 fprintf(fid,"%% %s : auto-generated file\n\n", OUTPUT_FILENAME);
235 fprintf(fid,"clear all\n");
236 fprintf(fid,"close all\n");
237
238 fprintf(fid,"k = %u;\n", k);
239 fprintf(fid,"m = %u;\n", m);
240 fprintf(fid,"num_symbols = %u;\n", num_symbols);
241 fprintf(fid,"num_samples = num_symbols*k;\n");
242
243 // save transmit matched-filter response
244 fprintf(fid,"hm_len = 2*k*m+1;\n");
245 fprintf(fid,"hm = zeros(1,hm_len);\n");
246 for (i=0; i<hm_len; i++)
247 fprintf(fid,"hm(%4u) = %12.4e;\n", i+1, hm[i]);
248
249 // save channel impulse response
250 fprintf(fid,"hc_len = %u;\n", hc_len);
251 fprintf(fid,"hc = zeros(1,hc_len);\n");
252 for (i=0; i<hc_len; i++)
253 fprintf(fid,"hc(%4u) = %12.4e + j*%12.4e;\n", i+1, crealf(hc[i]), cimagf(hc[i]));
254
255 // save equalizer response
256 fprintf(fid,"hp_len = %u;\n", hp_len);
257 fprintf(fid,"hp = zeros(1,hp_len);\n");
258 for (i=0; i<hp_len; i++)
259 fprintf(fid,"hp(%4u) = %12.4e + j*%12.4e;\n", i+1, crealf(hp[i]), cimagf(hp[i]));
260
261 // save sample sets
262 fprintf(fid,"x = zeros(1,num_samples);\n");
263 fprintf(fid,"y = zeros(1,num_samples);\n");
264 fprintf(fid,"z = zeros(1,num_samples);\n");
265 for (i=0; i<num_samples; i++) {
266 fprintf(fid,"x(%4u) = %12.4e + j*%12.4e;\n", i+1, crealf(x[i]), cimagf(x[i]));
267 fprintf(fid,"y(%4u) = %12.4e + j*%12.4e;\n", i+1, crealf(y[i]), cimagf(y[i]));
268 fprintf(fid,"z(%4u) = %12.4e + j*%12.4e;\n", i+1, crealf(z[i]), cimagf(z[i]));
269 }
270 fprintf(fid,"syms_rx = zeros(1,num_symbols);\n");
271 for (i=0; i<num_symbols; i++) {
272 fprintf(fid,"syms_rx(%4u) = %12.4e + j*%12.4e;\n", i+1, crealf(syms_rx[i]), cimagf(syms_rx[i]));
273 }
274
275 // plot time response
276 fprintf(fid,"t = 0:(num_samples-1);\n");
277 fprintf(fid,"tsym = 1:k:num_samples;\n");
278 fprintf(fid,"figure;\n");
279 fprintf(fid,"plot(t,real(z),...\n");
280 fprintf(fid," t(tsym),real(z(tsym)),'x');\n");
281
282 // plot constellation
283 fprintf(fid,"syms_rx_0 = syms_rx(1:(length(syms_rx)/2));\n");
284 fprintf(fid,"syms_rx_1 = syms_rx((length(syms_rx)/2):end);\n");
285 fprintf(fid,"figure;\n");
286 fprintf(fid,"plot(real(syms_rx_0),imag(syms_rx_0),'x','Color',[1 1 1]*0.7,...\n");
287 fprintf(fid," real(syms_rx_1),imag(syms_rx_1),'x','Color',[1 1 1]*0.0);\n");
288 fprintf(fid,"xlabel('In-Phase');\n");
289 fprintf(fid,"ylabel('Quadrature');\n");
290 fprintf(fid,"axis([-1 1 -1 1]*1.5);\n");
291 fprintf(fid,"axis square;\n");
292 fprintf(fid,"grid on;\n");
293
294 // compute composite response
295 fprintf(fid,"g = real(conv(conv(hm,hc),hp));\n");
296
297 // plot responses
298 fprintf(fid,"nfft = 1024;\n");
299 fprintf(fid,"f = [0:(nfft-1)]/nfft - 0.5;\n");
300 fprintf(fid,"Hm = 20*log10(abs(fftshift(fft(hm/k,nfft))));\n");
301 fprintf(fid,"Hc = 20*log10(abs(fftshift(fft(hc, nfft))));\n");
302 fprintf(fid,"Hp = 20*log10(abs(fftshift(fft(hp, nfft))));\n");
303 fprintf(fid,"G = 20*log10(abs(fftshift(fft(g/k, nfft))));\n");
304
305 fprintf(fid,"figure;\n");
306 fprintf(fid,"plot(f,Hm, f,Hc, f,Hp, f,G,'-k','LineWidth',2, [-0.5/k 0.5/k],[-6.026 -6.026],'or');\n");
307 fprintf(fid,"xlabel('Normalized Frequency');\n");
308 fprintf(fid,"ylabel('Power Spectral Density');\n");
309 fprintf(fid,"legend('transmit','channel','equalizer','composite','half-power points','location','northeast');\n");
310 fprintf(fid,"axis([-0.5 0.5 -12 8]);\n");
311 fprintf(fid,"grid on;\n");
312
313 fclose(fid);
314 printf("results written to '%s'\n", OUTPUT_FILENAME);
315
316 return 0;
317 }
318