1 /*
2 * Copyright (c) 2007 - 2017 Joseph Gaeddert
3 *
4 * Permission is hereby granted, free of charge, to any person obtaining a copy
5 * of this software and associated documentation files (the "Software"), to deal
6 * in the Software without restriction, including without limitation the rights
7 * to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
8 * copies of the Software, and to permit persons to whom the Software is
9 * furnished to do so, subject to the following conditions:
10 *
11 * The above copyright notice and this permission notice shall be included in
12 * all copies or substantial portions of the Software.
13 *
14 * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
15 * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
16 * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
17 * AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
18 * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
19 * OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
20 * THE SOFTWARE.
21 */
22
23 //
24 // gmskdem.c : Gauss minimum-shift keying modem
25 //
26
27 #include <stdio.h>
28 #include <assert.h>
29 #include <math.h>
30 #include <stdlib.h>
31
32 #include "liquid.internal.h"
33
34 #define DEBUG_GMSKDEM 0
35 #define DEBUG_GMSKDEM_FILENAME "gmskdem_internal_debug.m"
36 #define DEBUG_BUFFER_LEN (1000)
37
38 #define GMSKDEM_USE_EQUALIZER 0
39
40 void gmskdem_debug_print(gmskdem _q,
41 const char * _filename);
42
43 struct gmskdem_s {
44 unsigned int k; // samples/symbol
45 unsigned int m; // symbol delay
46 float BT; // bandwidth/time product
47 unsigned int h_len; // filter length
48 float * h; // pulse shaping filter
49
50 // filter object
51 #if GMSKDEM_USE_EQUALIZER
52 eqlms_rrrf eq; // receiver matched filter/equalizer
53 float k_inv; // 1 / k
54 #else
55 firfilt_rrrf filter; // receiver matched filter
56 #endif
57
58 float complex x_prime; // received signal state
59
60 // demodulated symbols counter
61 unsigned int num_symbols_demod;
62
63 #if DEBUG_GMSKDEM
64 windowf debug_mfout;
65 #endif
66 };
67
68 // create gmskdem object
69 // _k : samples/symbol
70 // _m : filter delay (symbols)
71 // _BT : excess bandwidth factor
gmskdem_create(unsigned int _k,unsigned int _m,float _BT)72 gmskdem gmskdem_create(unsigned int _k,
73 unsigned int _m,
74 float _BT)
75 {
76 if (_k < 2) {
77 fprintf(stderr,"error: gmskdem_create(), samples/symbol must be at least 2\n");
78 exit(1);
79 } else if (_m < 1) {
80 fprintf(stderr,"error: gmskdem_create(), symbol delay must be at least 1\n");
81 exit(1);
82 } else if (_BT <= 0.0f || _BT >= 1.0f) {
83 fprintf(stderr,"error: gmskdem_create(), bandwidth/time product must be in (0,1)\n");
84 exit(1);
85 }
86
87 // allocate memory for main object
88 gmskdem q = (gmskdem)malloc(sizeof(struct gmskdem_s));
89
90 // set properties
91 q->k = _k;
92 q->m = _m;
93 q->BT = _BT;
94
95 // allocate memory for filter taps
96 q->h_len = 2*(q->k)*(q->m)+1;
97 q->h = (float*) malloc(q->h_len * sizeof(float));
98
99 // compute filter coefficients
100 liquid_firdes_gmskrx(q->k, q->m, q->BT, 0.0f, q->h);
101
102 #if GMSKDEM_USE_EQUALIZER
103 // receiver matched filter/equalizer
104 q->eq = eqlms_rrrf_create_rnyquist(LIQUID_FIRFILT_GMSKRX,
105 q->k,
106 q->m,
107 q->BT,
108 0.0f);
109 eqlms_rrrf_set_bw(q->eq, 0.01f);
110 q->k_inv = 1.0f / (float)(q->k);
111 #else
112 // create filter object
113 q->filter = firfilt_rrrf_create(q->h, q->h_len);
114 #endif
115
116 // reset modem state
117 gmskdem_reset(q);
118
119 #if DEBUG_GMSKDEM
120 q->debug_mfout = windowf_create(DEBUG_BUFFER_LEN);
121 #endif
122
123 // return modem object
124 return q;
125 }
126
gmskdem_destroy(gmskdem _q)127 void gmskdem_destroy(gmskdem _q)
128 {
129 #if DEBUG_GMSKDEM
130 // print to external file
131 gmskdem_debug_print(_q, DEBUG_GMSKDEM_FILENAME);
132
133 // destroy debugging objects
134 windowf_destroy(_q->debug_mfout);
135 #endif
136
137 // destroy filter object
138 #if GMSKDEM_USE_EQUALIZER
139 eqlms_rrrf_destroy(_q->eq);
140 #else
141 firfilt_rrrf_destroy(_q->filter);
142 #endif
143
144 // free filter array
145 free(_q->h);
146
147 // free main object memory
148 free(_q);
149 }
150
gmskdem_print(gmskdem _q)151 void gmskdem_print(gmskdem _q)
152 {
153 printf("gmskdem [k=%u, m=%u, BT=%8.3f]\n", _q->k, _q->m, _q->BT);
154 #if GMSKDEM_USE_EQUALIZER
155 printf(" equalizer bandwidth : %12.8f\n", eqlms_rrrf_get_bw(_q->eq));
156 #endif
157 unsigned int i;
158 for (i=0; i<_q->h_len; i++)
159 printf(" hr(%4u) = %12.8f;\n", i+1, _q->h[i]);
160 }
161
gmskdem_reset(gmskdem _q)162 void gmskdem_reset(gmskdem _q)
163 {
164 // reset phase state
165 _q->x_prime = 0.0f;
166
167 // set demod. counter to zero
168 _q->num_symbols_demod = 0;
169
170 // clear filter buffer
171 #if GMSKDEM_USE_EQUALIZER
172 eqlms_rrrf_reset(_q->eq);
173 #else
174 firfilt_rrrf_reset(_q->filter);
175 #endif
176 }
177
178 // set equalizer bandwidth
gmskdem_set_eq_bw(gmskdem _q,float _bw)179 void gmskdem_set_eq_bw(gmskdem _q,
180 float _bw)
181 {
182 // validate input
183 if (_bw < 0.0f || _bw > 0.5f) {
184 fprintf(stderr,"error: gmskdem_set_eq_bw(), bandwith must be in [0,0.5]\n");
185 exit(1);
186 }
187
188 #if GMSKDEM_USE_EQUALIZER
189 // set internal equalizer bandwidth
190 eqlms_rrrf_set_bw(_q->eq, _bw);
191 #else
192 fprintf(stderr,"warning: gmskdem_set_eq_bw(), equalizer is disabled\n");
193 #endif
194 }
195
gmskdem_demodulate(gmskdem _q,float complex * _x,unsigned int * _s)196 void gmskdem_demodulate(gmskdem _q,
197 float complex * _x,
198 unsigned int * _s)
199 {
200 // increment symbol counter
201 _q->num_symbols_demod++;
202
203 // run matched filter
204 unsigned int i;
205 float phi;
206 float d_hat;
207 for (i=0; i<_q->k; i++) {
208 // compute phase difference
209 phi = cargf( conjf(_q->x_prime)*_x[i] );
210 _q->x_prime = _x[i];
211
212 // run through matched filter
213 #if GMSKDEM_USE_EQUALIZER
214 eqlms_rrrf_push(_q->eq, phi);
215 #else
216 firfilt_rrrf_push(_q->filter, phi);
217 #endif
218
219 #if DEBUG_GMSKDEM
220 // compute output
221 float d_tmp;
222 # if GMSKDEM_USE_EQUALIZER
223 eqlms_rrrf_execute(_q->eq, &d_tmp);
224 # else
225 firfilt_rrrf_execute(_q->filter, &d_tmp);
226 # endif
227 windowf_push(_q->debug_mfout, d_tmp);
228 #endif
229
230 // decimate by k
231 if ( i != 0 ) continue;
232
233 #if GMSKDEM_USE_EQUALIZER
234 // compute filter/equalizer output
235 eqlms_rrrf_execute(_q->eq, &d_hat);
236 #else
237 // compute filter output
238 firfilt_rrrf_execute(_q->filter, &d_hat);
239 #endif
240 }
241
242 // make decision
243 *_s = d_hat > 0.0f ? 1 : 0;
244
245 #if GMSKDEM_USE_EQUALIZER
246 // update equalizer, after appropriate delay
247 if (_q->num_symbols_demod >= 2*_q->m) {
248 // compute expected output, scaling by samples/symbol
249 float d_prime = d_hat > 0 ? _q->k_inv : -_q->k_inv;
250 eqlms_rrrf_step(_q->eq, d_prime, d_hat);
251 }
252 #endif
253 }
254
255 //
256 // output debugging file
257 //
gmskdem_debug_print(gmskdem _q,const char * _filename)258 void gmskdem_debug_print(gmskdem _q,
259 const char * _filename)
260 {
261 // open output filen for writing
262 FILE * fid = fopen(_filename,"w");
263 if (!fid) {
264 fprintf(stderr,"error: gmskdem_debug_print(), could not open '%s' for writing\n", _filename);
265 exit(1);
266 }
267 fprintf(fid,"%% %s : auto-generated file\n", _filename);
268 fprintf(fid,"clear all\n");
269 fprintf(fid,"close all\n");
270
271 #if DEBUG_GMSKDEM
272 //
273 unsigned int i;
274 float * r;
275 fprintf(fid,"n = %u;\n", DEBUG_BUFFER_LEN);
276 fprintf(fid,"k = %u;\n", _q->k);
277 fprintf(fid,"m = %u;\n", _q->m);
278 fprintf(fid,"t = [0:(n-1)]/k;\n");
279
280 // plot receive filter response
281 fprintf(fid,"ht = zeros(1,2*k*m+1);\n");
282 float ht[_q->h_len];
283 liquid_firdes_gmsktx(_q->k, _q->m, _q->BT, 0.0f, ht);
284 for (i=0; i<_q->h_len; i++)
285 fprintf(fid,"ht(%4u) = %12.4e;\n", i+1, ht[i]);
286 #if GMSKDEM_USE_EQUALIZER
287 float hr[_q->h_len];
288 eqlms_rrrf_get_weights(_q->eq, hr);
289 for (i=0; i<_q->h_len; i++)
290 fprintf(fid,"hr(%4u) = %12.4e * %u;\n", i+1, hr[i], _q->k);
291 #else
292 for (i=0; i<_q->h_len; i++)
293 fprintf(fid,"hr(%4u) = %12.4e;\n", i+1, _q->h[i]);
294 #endif
295 fprintf(fid,"hc = conv(ht,hr)/k;\n");
296 fprintf(fid,"nfft = 1024;\n");
297 fprintf(fid,"f = [0:(nfft-1)]/nfft - 0.5;\n");
298 fprintf(fid,"Ht = 20*log10(abs(fftshift(fft(ht/k, nfft))));\n");
299 fprintf(fid,"Hr = 20*log10(abs(fftshift(fft(hr/k, nfft))));\n");
300 fprintf(fid,"Hc = 20*log10(abs(fftshift(fft(hc/k, nfft))));\n");
301 fprintf(fid,"figure;\n");
302 fprintf(fid,"plot(f,Ht, f,Hr, f,Hc,'-k','LineWidth',2);\n");
303 fprintf(fid,"axis([-0.5 0.5 -50 10]);\n");
304 fprintf(fid,"xlabel('Normalized Frequency');\n");
305 fprintf(fid,"ylabel('Power Spectral Density [dB]');\n");
306 fprintf(fid,"legend('transmit','receive','composite',1);\n");
307 fprintf(fid,"grid on;\n");
308
309 fprintf(fid,"mfout = zeros(1,n);\n");
310 windowf_read(_q->debug_mfout, &r);
311 for (i=0; i<DEBUG_BUFFER_LEN; i++)
312 fprintf(fid,"mfout(%5u) = %12.4e;\n", i+1, r[i]);
313 fprintf(fid,"i0 = 1; %%mod(k+n,k)+k;\n");
314 fprintf(fid,"isym = i0:k:n;\n");
315 fprintf(fid,"figure;\n");
316 fprintf(fid,"plot(t,mfout,'-', t(isym),mfout(isym),'o','MarkerSize',4);\n");
317 fprintf(fid,"grid on;\n");
318 #endif
319
320 fclose(fid);
321 printf("gmskdem: internal debugging written to '%s'\n", _filename);
322 }
323