1 /*
2 * Copyright (c) 2007 - 2015 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 // Design root-Nyquist Kaiser filter
25 //
26 // References
27 // [Vaidyanathan:1993] Vaidyanathan, P. P., "Multirate Systems and
28 // Filter Banks," 1993, Prentice Hall, Section 3.2.1
29 //
30
31 #include <math.h>
32 #include <stdio.h>
33 #include <stdlib.h>
34
35 #include "liquid.internal.h"
36
37 #define DEBUG_RKAISER 0
38
39 // Design frequency-shifted root-Nyquist filter based on
40 // the Kaiser-windowed sinc.
41 //
42 // _k : filter over-sampling rate (samples/symbol)
43 // _m : filter delay (symbols)
44 // _beta : filter excess bandwidth factor (0,1)
45 // _dt : filter fractional sample delay
46 // _h : resulting filter [size: 2*_k*_m+1]
liquid_firdes_rkaiser(unsigned int _k,unsigned int _m,float _beta,float _dt,float * _h)47 void liquid_firdes_rkaiser(unsigned int _k,
48 unsigned int _m,
49 float _beta,
50 float _dt,
51 float * _h)
52 {
53 // validate input
54 if (_k < 2) {
55 fprintf(stderr,"error: liquid_firdes_rkaiser(), k must be at least 2\n");
56 exit(1);
57 } else if (_m < 1) {
58 fprintf(stderr,"error: liquid_firdes_rkaiser(), m must be at least 1\n");
59 exit(1);
60 } else if (_beta <= 0.0f || _beta >= 1.0f) {
61 fprintf(stderr,"error: liquid_firdes_rkaiser(), beta must be in (0,1)\n");
62 exit(1);
63 } else if (_dt < -1.0f || _dt > 1.0f) {
64 fprintf(stderr,"error: liquid_firdes_rkaiser(), dt must be in [-1,1]\n");
65 exit(1);
66 }
67
68 // simply call internal method and ignore output rho value
69 float rho;
70 //liquid_firdes_rkaiser_bisection(_k,_m,_beta,_dt,_h,&rho);
71 liquid_firdes_rkaiser_quadratic(_k,_m,_beta,_dt,_h,&rho);
72 }
73
74 // Design frequency-shifted root-Nyquist filter based on
75 // the Kaiser-windowed sinc using approximation for rho.
76 //
77 // _k : filter over-sampling rate (samples/symbol)
78 // _m : filter delay (symbols)
79 // _beta : filter excess bandwidth factor (0,1)
80 // _dt : filter fractional sample delay
81 // _h : resulting filter [size: 2*_k*_m+1]
liquid_firdes_arkaiser(unsigned int _k,unsigned int _m,float _beta,float _dt,float * _h)82 void liquid_firdes_arkaiser(unsigned int _k,
83 unsigned int _m,
84 float _beta,
85 float _dt,
86 float * _h)
87 {
88 // validate input
89 if (_k < 2) {
90 fprintf(stderr,"error: liquid_firdes_arkaiser(), k must be at least 2\n");
91 exit(1);
92 } else if (_m < 1) {
93 fprintf(stderr,"error: liquid_firdes_arkaiser(), m must be at least 1\n");
94 exit(1);
95 } else if (_beta <= 0.0f || _beta >= 1.0f) {
96 fprintf(stderr,"error: liquid_firdes_arkaiser(), beta must be in (0,1)\n");
97 exit(1);
98 } else if (_dt < -1.0f || _dt > 1.0f) {
99 fprintf(stderr,"error: liquid_firdes_arkaiser(), dt must be in [-1,1]\n");
100 exit(1);
101 }
102
103 #if 0
104 // compute bandwidth adjustment estimate
105 float rho_hat = rkaiser_approximate_rho(_m,_beta); // bandwidth correction factor
106 #else
107 // rho ~ c0 + c1*log(_beta) + c2*log^2(_beta)
108
109 // c0 ~ 0.762886 + 0.067663*log(m)
110 // c1 ~ 0.065515
111 // c2 ~ log( 1 - 0.088*m^-1.6 )
112
113 float c0 = 0.762886 + 0.067663*logf(_m);
114 float c1 = 0.065515;
115 float c2 = logf( 1 - 0.088*powf(_m,-1.6 ) );
116
117 float log_beta = logf(_beta);
118
119 float rho_hat = c0 + c1*log_beta + c2*log_beta*log_beta;
120
121 // ensure range is valid
122 if (rho_hat <= 0.0f || rho_hat >= 1.0f)
123 rho_hat = rkaiser_approximate_rho(_m,_beta);
124 #endif
125
126 unsigned int n=2*_k*_m+1; // filter length
127 float kf = (float)_k; // samples/symbol (float)
128 float del = _beta*rho_hat / kf; // transition bandwidth
129 float As = estimate_req_filter_As(del, n); // stop-band suppression
130 float fc = 0.5f*(1 + _beta*(1.0f-rho_hat))/kf; // filter cutoff
131
132 #if DEBUG_RKAISER
133 printf("rho-hat : %12.8f (compare to %12.8f)\n", rho_hat, rkaiser_approximate_rho(_m,_beta));
134 printf("fc : %12.8f\n", fc);
135 printf("delta-f : %12.8f\n", del);
136 printf("As : %12.8f dB\n", As);
137 printf("alpha : %12.8f\n", kaiser_beta_As(As));
138 #endif
139
140 // compute filter coefficients
141 liquid_firdes_kaiser(n,fc,As,_dt,_h);
142
143 // normalize coefficients
144 float e2 = 0.0f;
145 unsigned int i;
146 for (i=0; i<n; i++) e2 += _h[i]*_h[i];
147 for (i=0; i<n; i++) _h[i] *= sqrtf(_k/e2);
148 }
149
150 // Find approximate bandwidth adjustment factor rho based on
151 // filter delay and desired excess bandwdith factor.
152 //
153 // _m : filter delay (symbols)
154 // _beta : filter excess bandwidth factor (0,1)
rkaiser_approximate_rho(unsigned int _m,float _beta)155 float rkaiser_approximate_rho(unsigned int _m,
156 float _beta)
157 {
158 if ( _m < 1 ) {
159 fprintf(stderr,"error: rkaiser_approximate_rho(): m must be greater than 0\n");
160 exit(1);
161 } else if ( (_beta < 0.0f) || (_beta > 1.0f) ) {
162 fprintf(stderr,"error: rkaiser_approximate_rho(): beta must be in [0,1]\n");
163 exit(1);
164 } else;
165
166 // compute bandwidth adjustment estimate
167 float c0=0.0f, c1=0.0f, c2=0.0f;
168 switch (_m) {
169 case 1: c0=0.75749731; c1=0.06134303; c2=-0.08729663; break;
170 case 2: c0=0.81151861; c1=0.07437658; c2=-0.01427088; break;
171 case 3: c0=0.84249538; c1=0.07684185; c2=-0.00536879; break;
172 case 4: c0=0.86140782; c1=0.07144126; c2=-0.00558652; break;
173 case 5: c0=0.87457740; c1=0.06578694; c2=-0.00650447; break;
174 case 6: c0=0.88438797; c1=0.06074265; c2=-0.00736405; break;
175 case 7: c0=0.89216620; c1=0.05669236; c2=-0.00791222; break;
176 case 8: c0=0.89874983; c1=0.05361696; c2=-0.00815301; break;
177 case 9: c0=0.90460032; c1=0.05167952; c2=-0.00807893; break;
178 case 10: c0=0.91034430; c1=0.05130753; c2=-0.00746192; break;
179 case 11: c0=0.91587675; c1=0.05180436; c2=-0.00670711; break;
180 case 12: c0=0.92121875; c1=0.05273801; c2=-0.00588351; break;
181 case 13: c0=0.92638195; c1=0.05400764; c2=-0.00508452; break;
182 case 14: c0=0.93123555; c1=0.05516163; c2=-0.00437306; break;
183 case 15: c0=0.93564993; c1=0.05596561; c2=-0.00388152; break;
184 case 16: c0=0.93976742; c1=0.05662274; c2=-0.00348280; break;
185 case 17: c0=0.94351703; c1=0.05694120; c2=-0.00318821; break;
186 case 18: c0=0.94557273; c1=0.05227591; c2=-0.00400676; break;
187 case 19: c0=0.95001614; c1=0.05681641; c2=-0.00300628; break;
188 case 20: c0=0.95281708; c1=0.05637607; c2=-0.00304790; break;
189 case 21: c0=0.95536256; c1=0.05575880; c2=-0.00312988; break;
190 case 22: c0=0.95754206; c1=0.05426060; c2=-0.00385945; break;
191 default:
192 c0 = 0.056873*logf(_m+1e-3f) + 0.781388;
193 c1 = 0.05426f;
194 c2 = -0.00386f;
195 }
196
197 float b = logf(_beta);
198 float rho_hat = c0 + c1*b + c2*b*b;
199
200 // ensure estimate is in [0,1]
201 if (rho_hat < 0.0f) {
202 rho_hat = 0.0f;
203 } else if (rho_hat > 1.0f) {
204 rho_hat = 1.0f;
205 }
206
207 return rho_hat;
208 }
209
210 // Design frequency-shifted root-Nyquist filter based on
211 // the Kaiser-windowed sinc.
212 //
213 // _k : filter over-sampling rate (samples/symbol)
214 // _m : filter delay (symbols)
215 // _beta : filter excess bandwidth factor (0,1)
216 // _dt : filter fractional sample delay
217 // _h : resulting filter [size: 2*_k*_m+1]
218 // _rho : transition bandwidth adjustment, 0 < _rho < 1
liquid_firdes_rkaiser_bisection(unsigned int _k,unsigned int _m,float _beta,float _dt,float * _h,float * _rho)219 void liquid_firdes_rkaiser_bisection(unsigned int _k,
220 unsigned int _m,
221 float _beta,
222 float _dt,
223 float * _h,
224 float * _rho)
225 {
226 if ( _k < 1 ) {
227 fprintf(stderr,"error: liquid_firdes_rkaiser_bisection(): k must be greater than 0\n");
228 exit(1);
229 } else if ( _m < 1 ) {
230 fprintf(stderr,"error: liquid_firdes_rkaiser_bisection(): m must be greater than 0\n");
231 exit(1);
232 } else if ( (_beta < 0.0f) || (_beta > 1.0f) ) {
233 fprintf(stderr,"error: liquid_firdes_rkaiser_bisection(): beta must be in [0,1]\n");
234 exit(1);
235 } else;
236
237 // algorithm:
238 // 1. choose three initial points [x0, x1, x2] where x0 < x1 < x2
239 // 2. choose xa = 0.5*(x0 + x1), bisection of [x0 x1]
240 // 3. choose xb = 0.5*(x1 + x2), bisection of [x1 x2]
241 // 4. x0 < xa < x1 < xb < x2
242 // 5. evaluate points to obtain [y0, ya, y1, yb, y2]
243 // 6. find minimum of three midpoints [ya, y1, yb]
244 // 7. update initial points [x0, x1, x2] and go to step 2
245
246 unsigned int i;
247
248 unsigned int n=2*_k*_m+1; // filter length
249
250 // compute bandwidth adjustment estimate
251 float rho_hat = rkaiser_approximate_rho(_m,_beta);
252
253 // bandwidth adjustment
254 float x0 = 0.5f * rho_hat; // lower bound, old: rho_hat*0.9f;
255 float x1 = rho_hat; // midpoint: use initial estimate
256 float x2 = 1.0f; // upper bound, old: 1.0f - 0.9f*(1.0f-rho_hat);
257 //x1 = 0.5f*(x0 + x2); // bisect [x0,x1]
258
259 // evaluate performance (ISI) of each bandwidth adjustment
260 float y0 = liquid_firdes_rkaiser_internal_isi(_k,_m,_beta,_dt,x0,_h);
261 float y1 = liquid_firdes_rkaiser_internal_isi(_k,_m,_beta,_dt,x1,_h);
262 float y2 = liquid_firdes_rkaiser_internal_isi(_k,_m,_beta,_dt,x2,_h);
263
264 // run parabolic search to find bandwidth adjustment x_hat which
265 // minimizes the inter-symbol interference of the filter
266 unsigned int p, pmax=14;
267 float x_hat = rho_hat;
268 float xa, xb;
269 float ya, yb;
270 #if DEBUG_RKAISER
271 FILE * fid = fopen("rkaiser_debug.m", "w");
272 fprintf(fid,"clear all;\n");
273 fprintf(fid,"close all;\n");
274 fprintf(fid,"x = [%12.4e %12.4e %12.4e];\n", x0, x1, x2);
275 fprintf(fid,"y = [%12.4e %12.4e %12.4e];\n", y0, y1, y2);
276 #endif
277 for (p=0; p<pmax; p++) {
278 // check bounding conditions: y1 should be less than y0 and y2
279 if (y1 > y0 || y1 > y2)
280 fprintf(stderr,"warning: liquid_firdes_rkaiser_bisection(): bounding region is ill-conditioned\n");
281
282 // choose midway points xa, xb and compute ISI
283 xa = 0.5f*(x0 + x1); // bisect [x0,x1]
284 xb = 0.5f*(x1 + x2); // bisect [x1,x2]
285 ya = liquid_firdes_rkaiser_internal_isi(_k,_m,_beta,_dt,xa,_h);
286 yb = liquid_firdes_rkaiser_internal_isi(_k,_m,_beta,_dt,xb,_h);
287
288 #if DEBUG_RKAISER
289 fprintf(fid,"x = [x %12.4e %12.4e];\n", xa, xb);
290 fprintf(fid,"y = [y %12.4e %12.4e];\n", ya, yb);
291 #endif
292
293 // find the minimum of [ya,y1,yb], update bounds
294 if (y1 < ya && y1 < yb) {
295 x0 = xa; y0 = ya;
296 x2 = xb; y2 = yb;
297 } else if (ya < yb) {
298 x2 = x1; y2 = y1;
299 x1 = xa; y1 = ya;
300 } else {
301 x0 = x1; y0 = y1;
302 x1 = xb; y1 = yb;
303 }
304
305 x_hat = x1;
306 #if DEBUG_RKAISER
307 float y_hat = liquid_firdes_rkaiser_internal_isi(_k,_m,_beta,_dt,x_hat,_h);
308 printf(" %4u : rho=%12.8f, isi=%12.6f dB\n", p+1, x_hat, 20*log10f(y_hat));
309 #endif
310 };
311 #if DEBUG_RKAISER
312 fprintf(fid,"figure;\n");
313 fprintf(fid,"plot(x,20*log10(y),'x');\n");
314 fclose(fid);
315 printf("results written to rkaiser_debug.m\n");
316 #endif
317
318 // re-design filter with optimal value for rho
319 liquid_firdes_rkaiser_internal_isi(_k,_m,_beta,_dt,x_hat,_h);
320
321 // normalize filter magnitude
322 float e2 = 0.0f;
323 for (i=0; i<n; i++) e2 += _h[i]*_h[i];
324 for (i=0; i<n; i++) _h[i] *= sqrtf(_k/e2);
325
326 // save trasition bandwidth adjustment
327 *_rho = x_hat;
328 }
329
330 // Design frequency-shifted root-Nyquist filter based on
331 // the Kaiser-windowed sinc using the quadratic search method.
332 //
333 // _k : filter over-sampling rate (samples/symbol)
334 // _m : filter delay (symbols)
335 // _beta : filter excess bandwidth factor (0,1)
336 // _dt : filter fractional sample delay
337 // _h : resulting filter [size: 2*_k*_m+1]
338 // _rho : transition bandwidth adjustment, 0 < _rho < 1
liquid_firdes_rkaiser_quadratic(unsigned int _k,unsigned int _m,float _beta,float _dt,float * _h,float * _rho)339 void liquid_firdes_rkaiser_quadratic(unsigned int _k,
340 unsigned int _m,
341 float _beta,
342 float _dt,
343 float * _h,
344 float * _rho)
345 {
346 if ( _k < 1 ) {
347 fprintf(stderr,"error: liquid_firdes_rkaiser_quadratic(): k must be greater than 0\n");
348 exit(1);
349 } else if ( _m < 1 ) {
350 fprintf(stderr,"error: liquid_firdes_rkaiser_quadratic(): m must be greater than 0\n");
351 exit(1);
352 } else if ( (_beta < 0.0f) || (_beta > 1.0f) ) {
353 fprintf(stderr,"error: liquid_firdes_rkaiser_quadratic(): beta must be in [0,1]\n");
354 exit(1);
355 } else;
356
357 // algorithm:
358 // 1. choose initial bounding points [x0,x2] where x0 < x2
359 // 2. choose x1 as bisection of [x0,x2]: x1 = 0.5*(x0+x2)
360 // 3. choose x_hat as solution to quadratic equation (x0,y0), (x1,y1), (x2,y2)
361 // 4. re-select boundary: (x0,y0) <- (x1,y1) if x_hat > x1
362 // (x2,y2) <- (x1,y1) otherwise
363 // 5. go to step 2
364
365 unsigned int i;
366
367 unsigned int n=2*_k*_m+1; // filter length
368
369 // compute bandwidth adjustment estimate
370 float rho_hat = rkaiser_approximate_rho(_m,_beta);
371 float rho_opt = rho_hat;
372
373 // bandwidth adjustment
374 float x0; // lower bound
375 float x1 = rho_hat; // initial estimate
376 float x2; // upper bound
377
378 // evaluate performance (ISI) of each bandwidth adjustment
379 float y0;
380 float y1;
381 float y2;
382 float y_opt = 0.0f;
383
384 // run parabolic search to find bandwidth adjustment x_hat which
385 // minimizes the inter-symbol interference of the filter
386 unsigned int p, pmax=14;
387 float x_hat = rho_hat;
388 float dx = 0.2f; // bounding size
389 float del = 0.0f; // computed step size
390 float tol = 1e-6f; // tolerance
391 #if DEBUG_RKAISER
392 FILE * fid = fopen("rkaiser_debug.m", "w");
393 fprintf(fid,"clear all;\n");
394 fprintf(fid,"close all;\n");
395 fprintf(fid,"x = [];\n");
396 fprintf(fid,"y = [];\n");
397 #endif
398 for (p=0; p<pmax; p++) {
399 // choose boundary points
400 x0 = x1 - dx;
401 x2 = x1 + dx;
402
403 // ensure boundary points are valid
404 if (x0 <= 0.0f) x0 = 0.01f;
405 if (x2 >= 1.0f) x2 = 0.99f;
406
407 // evaluate all points
408 y0 = liquid_firdes_rkaiser_internal_isi(_k,_m,_beta,_dt,x0,_h);
409 y1 = liquid_firdes_rkaiser_internal_isi(_k,_m,_beta,_dt,x1,_h);
410 y2 = liquid_firdes_rkaiser_internal_isi(_k,_m,_beta,_dt,x2,_h);
411
412 // save optimum
413 if (p==0 || y1 < y_opt) {
414 rho_opt = x1;
415 y_opt = y1;
416 }
417
418 #if DEBUG_RKAISER
419 fprintf(fid,"x = [x %12.4e %12.4e %12.4e];\n", x0, x1, x2);
420 fprintf(fid,"y = [y %12.4e %12.4e %12.4e];\n", y0, y1, y2);
421 //printf(" %4u : [%8.4f %8.4f %8.4f] rho=%12.8f, isi=%12.6f dB\n", p+1, x0, x1, x2, x1, 20*log10f(y1));
422 printf(" %4u : [%8.4f,%8.2f] [%8.4f,%8.2f] [%8.4f,%8.2f]\n",
423 p+1,
424 x0, 20*log10f(y0),
425 x1, 20*log10f(y1),
426 x2, 20*log10f(y2));
427 #endif
428 // compute minimum of quadratic function
429 double ta = y0*(x1*x1 - x2*x2) +
430 y1*(x2*x2 - x0*x0) +
431 y2*(x0*x0 - x1*x1);
432
433 double tb = y0*(x1 - x2) +
434 y1*(x2 - x0) +
435 y2*(x0 - x1);
436
437 // update estimate
438 x_hat = 0.5f * ta / tb;
439
440 // ensure x_hat is within boundary (this will fail if y1 > y0 || y1 > y2)
441 if (x_hat < x0 || x_hat > x2) {
442 //fprintf(stderr,"warning: liquid_firdes_rkaiser_quadratic(), quadratic minimum outside boundary\n");
443 break;
444 }
445
446 // break if step size is sufficiently small
447 del = x_hat - x1;
448 if (p > 3 && fabsf(del) < tol)
449 break;
450
451 // update estimate, reduce bound
452 x1 = x_hat;
453 dx *= 0.5f;
454 };
455 #if DEBUG_RKAISER
456 fprintf(fid,"figure;\n");
457 fprintf(fid,"plot(x,20*log10(y),'x');\n");
458 fclose(fid);
459 printf("results written to rkaiser_debug.m\n");
460 #endif
461
462 // re-design filter with optimal value for rho
463 liquid_firdes_rkaiser_internal_isi(_k,_m,_beta,_dt,rho_opt,_h);
464
465 // normalize filter magnitude
466 float e2 = 0.0f;
467 for (i=0; i<n; i++) e2 += _h[i]*_h[i];
468 for (i=0; i<n; i++) _h[i] *= sqrtf(_k/e2);
469
470 // save trasition bandwidth adjustment
471 *_rho = rho_opt;
472 }
473
474 // compute filter coefficients and determine resulting ISI
475 //
476 // _k : filter over-sampling rate (samples/symbol)
477 // _m : filter delay (symbols)
478 // _beta : filter excess bandwidth factor (0,1)
479 // _dt : filter fractional sample delay
480 // _rho : transition bandwidth adjustment, 0 < _rho < 1
481 // _h : filter buffer [size: 2*_k*_m+1]
liquid_firdes_rkaiser_internal_isi(unsigned int _k,unsigned int _m,float _beta,float _dt,float _rho,float * _h)482 float liquid_firdes_rkaiser_internal_isi(unsigned int _k,
483 unsigned int _m,
484 float _beta,
485 float _dt,
486 float _rho,
487 float * _h)
488 {
489 // validate input
490 if (_rho < 0.0f) {
491 fprintf(stderr,"warning: liquid_firdes_rkaiser_internal_isi(), rho < 0\n");
492 } else if (_rho > 1.0f) {
493 fprintf(stderr,"warning: liquid_firdes_rkaiser_internal_isi(), rho > 1\n");
494 }
495
496 unsigned int n=2*_k*_m+1; // filter length
497 float kf = (float)_k; // samples/symbol (float)
498 float del = _beta*_rho / kf; // transition bandwidth
499 float As = estimate_req_filter_As(del, n); // stop-band suppression
500 float fc = 0.5f*(1 + _beta*(1.0f-_rho))/kf; // filter cutoff
501
502 // evaluate performance (ISI)
503 float isi_max;
504 float isi_rms;
505
506 // compute filter
507 liquid_firdes_kaiser(n,fc,As,_dt,_h);
508
509 // compute filter ISI
510 liquid_filter_isi(_h,_k,_m,&isi_rms,&isi_max);
511
512 // return RMS of ISI
513 return isi_rms;
514 }
515
516
517