1 /*
2 * Copyright (c) 2007 - 2019 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 // Finite impulse response filter design
25 //
26 // References:
27 // [Herrmann:1973] O. Herrmann, L. R. Rabiner, and D. S. K. Chan,
28 // "Practical design rules for optimum finite impulse response
29 // lowpass digital filters," Bell Syst. Tech. Journal, vol. 52,
30 // pp. 769--99, July-Aug. 1973
31 // [Vaidyanathan:1993] Vaidyanathan, P. P., "Multirate Systems and
32 // Filter Banks," 1993, Prentice Hall, Section 3.2.1
33
34 #include <stdio.h>
35 #include <stdlib.h>
36 #include <string.h>
37 #include <math.h>
38
39 #include "liquid.internal.h"
40
41 // method to estimate required filter length
42 #define ESTIMATE_REQ_FILTER_LEN_METHOD_KAISER (0)
43 #define ESTIMATE_REQ_FILTER_LEN_METHOD_HERRMANN (1)
44
45 // select filter estimate method
46 #define ESTIMATE_REQ_FILTER_LEN_METHOD (0)
47
48 // esimate required filter length given transition bandwidth and
49 // stop-band attenuation
50 // _df : transition bandwidth (0 < _df < 0.5)
51 // _As : stopband suppression level [dB] (_As > 0)
estimate_req_filter_len(float _df,float _As)52 unsigned int estimate_req_filter_len(float _df,
53 float _As)
54 {
55 if (_df > 0.5f || _df <= 0.0f) {
56 fprintf(stderr,"error: estimate_req_filter_len(), invalid bandwidth : %f\n", _df);
57 exit(1);
58 } else if (_As <= 0.0f) {
59 fprintf(stderr,"error: estimate_req_filter_len(), invalid stopband level : %f\n", _As);
60 exit(1);
61 }
62
63 // compute filter length estimate
64 #if ESTIMATE_REQ_FILTER_LEN_METHOD == ESTIMATE_REQ_FILTER_LEN_METHOD_KAISER
65 // use Kaiser's estimate
66 unsigned int h_len = (unsigned int) estimate_req_filter_len_Kaiser(_df,_As);
67 #elif ESTIMATE_REQ_FILTER_LEN_METHOD == ESTIMATE_REQ_FILTER_LEN_METHOD_HERRMANN
68 // use Herrmann's estimate
69 unsigned int h_len = (unsigned int) estimate_req_filter_len_Herrmann(_df,_As);
70 #else
71 # error "invalid required filter length estimation method"
72 #endif
73
74 return h_len;
75 }
76
77 // estimate filter stop-band attenuation given
78 // _df : transition bandwidth (0 < _b < 0.5)
79 // _N : filter length
estimate_req_filter_As(float _df,unsigned int _N)80 float estimate_req_filter_As(float _df,
81 unsigned int _N)
82 {
83 // run search for stop-band attenuation which gives these results
84 float As0 = 0.01f; // lower bound
85 float As1 = 200.0f; // upper bound
86
87 float As_hat = 0.0f; // stop-band attenuation estimate
88 float N_hat = 0.0f; // filter length estimate
89
90 // perform simple bisection search
91 unsigned int num_iterations = 20;
92 unsigned int i;
93 for (i=0; i<num_iterations; i++) {
94 // bisect limits
95 As_hat = 0.5f*(As1 + As0);
96 #if ESTIMATE_REQ_FILTER_LEN_METHOD == ESTIMATE_REQ_FILTER_LEN_METHOD_KAISER
97 N_hat = estimate_req_filter_len_Kaiser(_df, As_hat);
98 #elif ESTIMATE_REQ_FILTER_LEN_METHOD == ESTIMATE_REQ_FILTER_LEN_METHOD_HERRMANN
99 N_hat = estimate_req_filter_len_Herrmann(_df, As_hat);
100 #else
101 # error "invalid required filter length estimation method"
102 #endif
103
104 //printf("range[%8.2f, %8.2f] As-hat=%8.2fdB, N=%8.2f (target: %3u taps)\n",
105 // As0, As1, As_hat, N_hat, _N);
106
107 // update limits
108 if (N_hat < (float)_N) {
109 As0 = As_hat;
110 } else {
111 As1 = As_hat;
112 }
113 }
114 return As_hat;
115 }
116
117 // estimate filter transition bandwidth given
118 // _As : stop-band attenuation [dB], _As > 0
119 // _N : filter length
estimate_req_filter_df(float _As,unsigned int _N)120 float estimate_req_filter_df(float _As,
121 unsigned int _N)
122 {
123 // run search for stop-band attenuation which gives these results
124 float df0 = 1e-3f; // lower bound
125 float df1 = 0.499f; // upper bound
126
127 float df_hat = 0.0f; // stop-band attenuation estimate
128 float N_hat = 0.0f; // filter length estimate
129
130 // perform simple bisection search
131 unsigned int num_iterations = 20;
132 unsigned int i;
133 for (i=0; i<num_iterations; i++) {
134 // bisect limits
135 df_hat = 0.5f*(df1 + df0);
136 #if ESTIMATE_REQ_FILTER_LEN_METHOD == ESTIMATE_REQ_FILTER_LEN_METHOD_KAISER
137 N_hat = estimate_req_filter_len_Kaiser(df_hat, _As);
138 #elif ESTIMATE_REQ_FILTER_LEN_METHOD == ESTIMATE_REQ_FILTER_LEN_METHOD_HERRMANN
139 N_hat = estimate_req_filter_len_Herrmann(df_hat, _As);
140 #else
141 # error "invalid required filter length estimation method"
142 #endif
143
144 //printf("range[%8.5f, %8.5f] df-hat=%8.5fdB, N=%8.2f (target: %3u taps)\n",
145 // df0, df1, df_hat, N_hat, _N);
146
147 // update limits
148 if (N_hat < (float)_N) {
149 df1 = df_hat;
150 } else {
151 df0 = df_hat;
152 }
153 }
154 return df_hat;
155
156 }
157
158
159 // esimate required filter length given transition bandwidth and
160 // stop-band attenuation (algorithm from [Vaidyanathan:1993])
161 // _df : transition bandwidth (0 < _df < 0.5)
162 // _As : stop-band attenuation [dB] (As > 0)
estimate_req_filter_len_Kaiser(float _df,float _As)163 float estimate_req_filter_len_Kaiser(float _df,
164 float _As)
165 {
166 if (_df > 0.5f || _df <= 0.0f) {
167 fprintf(stderr,"error: estimate_req_filter_len_Kaiser(), invalid bandwidth : %f\n", _df);
168 exit(1);
169 } else if (_As <= 0.0f) {
170 fprintf(stderr,"error: estimate_req_filter_len(), invalid stopband level : %f\n", _As);
171 exit(1);
172 }
173
174 // compute filter length estimate
175 return (_As - 7.95f)/(14.26f*_df);
176 }
177
178
179 // esimate required filter length given transition bandwidth and
180 // stop-band attenuation (algorithm from [Herrmann:1973])
181 // _df : transition bandwidth (0 < _df < 0.5)
182 // _As : stop-band attenuation [dB] (As > 0)
estimate_req_filter_len_Herrmann(float _df,float _As)183 float estimate_req_filter_len_Herrmann(float _df,
184 float _As)
185 {
186 if (_df > 0.5f || _df <= 0.0f) {
187 fprintf(stderr,"error: estimate_req_filter_len_Herrmann(), invalid bandwidth : %f\n", _df);
188 exit(1);
189 } else if (_As <= 0.0f) {
190 fprintf(stderr,"error: estimate_req_filter_len(), invalid stopband level : %f\n", _As);
191 exit(1);
192 }
193
194 // Gaeddert's revisions:
195 if (_As > 105.0f)
196 return estimate_req_filter_len_Kaiser(_df,_As);
197
198 _As += 7.4f;
199
200 // compute delta_1, delta_2
201 float d1, d2;
202 d1 = d2 = powf(10.0, -_As/20.0);
203
204 // compute log of delta_1, delta_2
205 float t1 = log10f(d1);
206 float t2 = log10f(d2);
207
208 // compute D_infinity(delta_1, delta_2)
209 float Dinf = (0.005309f*t1*t1 + 0.07114f*t1 - 0.4761f)*t2 -
210 (0.002660f*t1*t1 + 0.59410f*t1 + 0.4278f);
211
212 // compute f(delta_1, delta_2)
213 float f = 11.012f + 0.51244f*(t1-t2);
214
215 // compute filter length estimate
216 float h_len = (Dinf - f*_df*_df) / _df + 1.0f;
217
218 return h_len;
219 }
220
221 // returns the Kaiser window beta factor give the filter's target
222 // stop-band attenuation (As) [Vaidyanathan:1993]
223 // _As : target filter's stop-band attenuation [dB], _As > 0
kaiser_beta_As(float _As)224 float kaiser_beta_As(float _As)
225 {
226 _As = fabsf(_As);
227 float beta;
228 if (_As > 50.0f)
229 beta = 0.1102f*(_As - 8.7f);
230 else if (_As > 21.0f)
231 beta = 0.5842*powf(_As - 21, 0.4f) + 0.07886f*(_As - 21);
232 else
233 beta = 0.0f;
234
235 return beta;
236 }
237
238 // Design FIR using kaiser window
239 // _n : filter length, _n > 0
240 // _fc : cutoff frequency, 0 < _fc < 0.5
241 // _As : stop-band attenuation [dB], _As > 0
242 // _mu : fractional sample offset, -0.5 < _mu < 0.5
243 // _h : output coefficient buffer, [size: _n x 1]
liquid_firdes_kaiser(unsigned int _n,float _fc,float _As,float _mu,float * _h)244 void liquid_firdes_kaiser(unsigned int _n,
245 float _fc,
246 float _As,
247 float _mu,
248 float *_h)
249 {
250 // validate inputs
251 if (_mu < -0.5f || _mu > 0.5f) {
252 fprintf(stderr,"error: liquid_firdes_kaiser(), _mu (%12.4e) out of range [-0.5,0.5]\n", _mu);
253 exit(1);
254 } else if (_fc < 0.0f || _fc > 0.5f) {
255 fprintf(stderr,"error: liquid_firdes_kaiser(), cutoff frequency (%12.4e) out of range (0, 0.5)\n", _fc);
256 exit(1);
257 } else if (_n == 0) {
258 fprintf(stderr,"error: liquid_firdes_kaiser(), filter length must be greater than zero\n");
259 exit(1);
260 }
261
262 // choose kaiser beta parameter (approximate)
263 float beta = kaiser_beta_As(_As);
264
265 float t, h1, h2;
266 unsigned int i;
267 for (i=0; i<_n; i++) {
268 t = (float)i - (float)(_n-1)/2 + _mu;
269
270 // sinc prototype
271 h1 = sincf(2.0f*_fc*t);
272
273 // kaiser window
274 h2 = kaiser(i,_n,beta,_mu);
275
276 //printf("t = %f, h1 = %f, h2 = %f\n", t, h1, h2);
277
278 // composite
279 _h[i] = h1*h2;
280 }
281 }
282
283 // Design finite impulse response notch filter
284 // _m : filter semi-length, m in [1,1000]
285 // _f0 : filter notch frequency (normalized), -0.5 <= _fc <= 0.5
286 // _As : stop-band attenuation [dB], _As > 0
287 // _h : output coefficient buffer, [size: 2*_m+1 x 1]
liquid_firdes_notch(unsigned int _m,float _f0,float _As,float * _h)288 void liquid_firdes_notch(unsigned int _m,
289 float _f0,
290 float _As,
291 float * _h)
292 {
293 // validate inputs
294 if (_m < 1 || _m > 1000) {
295 fprintf(stderr,"error: liquid_firdes_notch(), _m (%12u) out of range [1,1000]\n", _m);
296 exit(1);
297 } else if (_f0 < -0.5f || _f0 > 0.5f) {
298 fprintf(stderr,"error: liquid_firdes_notch(), notch frequency (%12.4e) must be in [-0.5,0.5]\n", _f0);
299 exit(1);
300 } else if (_As <= 0.0f) {
301 fprintf(stderr,"error: liquid_firdes_notch(), stop-band suppression (%12.4e) must be greater than zero\n", _As);
302 exit(1);
303 }
304
305 // choose kaiser beta parameter (approximate)
306 float beta = kaiser_beta_As(_As);
307
308 // design filter
309 unsigned int h_len = 2*_m+1;
310 unsigned int i;
311 float scale = 0.0f;
312 for (i=0; i<h_len; i++) {
313 // tone at carrier frequency
314 float p = -cosf(2.0f*M_PI*_f0*((float)(i) - (float)_m));
315
316 // window
317 float w = kaiser(i,h_len,beta,0);
318
319 // save un-normalized filter
320 _h[i] = p*w;
321
322 // accumulate scale
323 scale += _h[i] * p;
324 }
325
326 // normalize
327 for (i=0; i<h_len; i++)
328 _h[i] /= scale;
329
330 // add impulse
331 _h[_m] += 1.0f;
332 }
333
334 // Design (root-)Nyquist filter from prototype
335 // _type : filter type (e.g. LIQUID_FIRFILT_RRRC)
336 // _k : samples/symbol
337 // _m : symbol delay
338 // _beta : excess bandwidth factor, _beta in [0,1]
339 // _dt : fractional sample delay
340 // _h : output coefficient buffer (length: 2*k*m+1)
liquid_firdes_prototype(liquid_firfilt_type _type,unsigned int _k,unsigned int _m,float _beta,float _dt,float * _h)341 void liquid_firdes_prototype(liquid_firfilt_type _type,
342 unsigned int _k,
343 unsigned int _m,
344 float _beta,
345 float _dt,
346 float * _h)
347 {
348 // compute filter parameters
349 unsigned int h_len = 2*_k*_m + 1; // length
350 float fc = 0.5f / (float)_k; // cut-off frequency
351 float df = _beta / (float)_k; // transition bandwidth
352 float As = estimate_req_filter_As(df,h_len); // stop-band attenuation
353
354 // Parks-McClellan algorithm parameters
355 float bands[6] = { 0.0f, fc-0.5f*df,
356 fc, fc,
357 fc+0.5f*df, 0.5f};
358 float des[3] = { (float)_k, 0.5f*_k, 0.0f };
359 float weights[3] = {1.0f, 1.0f, 1.0f};
360 liquid_firdespm_wtype wtype[3] = { LIQUID_FIRDESPM_FLATWEIGHT,
361 LIQUID_FIRDESPM_FLATWEIGHT,
362 LIQUID_FIRDESPM_FLATWEIGHT};
363
364 switch (_type) {
365
366 // Nyquist filter prototypes
367
368 case LIQUID_FIRFILT_KAISER:
369 liquid_firdes_kaiser(h_len, fc, As, _dt, _h);
370 break;
371 case LIQUID_FIRFILT_PM:
372 // WARNING: input timing offset is ignored here
373 firdespm_run(h_len, 3, bands, des, weights, wtype, LIQUID_FIRDESPM_BANDPASS, _h);
374 break;
375 case LIQUID_FIRFILT_RCOS:
376 liquid_firdes_rcos(_k, _m, _beta, _dt, _h);
377 break;
378 case LIQUID_FIRFILT_FEXP:
379 liquid_firdes_fexp(_k, _m, _beta, _dt, _h);
380 break;
381 case LIQUID_FIRFILT_FSECH:
382 liquid_firdes_fsech(_k, _m, _beta, _dt, _h);
383 break;
384 case LIQUID_FIRFILT_FARCSECH:
385 liquid_firdes_farcsech(_k, _m, _beta, _dt, _h);
386 break;
387
388 // root-Nyquist filter prototypes
389
390 case LIQUID_FIRFILT_ARKAISER:
391 liquid_firdes_arkaiser(_k, _m, _beta, _dt, _h);
392 break;
393 case LIQUID_FIRFILT_RKAISER:
394 liquid_firdes_rkaiser(_k, _m, _beta, _dt, _h);
395 break;
396 case LIQUID_FIRFILT_RRC:
397 liquid_firdes_rrcos(_k, _m, _beta, _dt, _h);
398 break;
399 case LIQUID_FIRFILT_hM3:
400 liquid_firdes_hM3(_k, _m, _beta, _dt, _h);
401 break;
402 case LIQUID_FIRFILT_GMSKTX:
403 liquid_firdes_gmsktx(_k, _m, _beta, _dt, _h);
404 break;
405 case LIQUID_FIRFILT_GMSKRX:
406 liquid_firdes_gmskrx(_k, _m, _beta, _dt, _h);
407 break;
408 case LIQUID_FIRFILT_RFEXP:
409 liquid_firdes_rfexp(_k, _m, _beta, _dt, _h);
410 break;
411 case LIQUID_FIRFILT_RFSECH:
412 liquid_firdes_rfsech(_k, _m, _beta, _dt, _h);
413 break;
414 case LIQUID_FIRFILT_RFARCSECH:
415 liquid_firdes_rfarcsech(_k, _m, _beta, _dt, _h);
416 break;
417 default:
418 fprintf(stderr,"error: liquid_firdes_prototype(), invalid root-Nyquist filter type '%d'\n", _type);
419 exit(1);
420 }
421 }
422
423
424 // Design FIR doppler filter
425 // _n : filter length
426 // _fd : normalized doppler frequency (0 < _fd < 0.5)
427 // _K : Rice fading factor (K >= 0)
428 // _theta : LoS component angle of arrival
429 // _h : output coefficient buffer
liquid_firdes_doppler(unsigned int _n,float _fd,float _K,float _theta,float * _h)430 void liquid_firdes_doppler(unsigned int _n,
431 float _fd,
432 float _K,
433 float _theta,
434 float * _h)
435 {
436 float t, J, r, w;
437 float beta = 4; // kaiser window parameter
438 unsigned int i;
439 for (i=0; i<_n; i++) {
440 // time sample
441 t = (float)i - (float)(_n-1)/2;
442
443 // Bessel
444 J = 1.5*liquid_besselj0f(fabsf((float)(2*M_PI*_fd*t)));
445
446 // Rice-K component
447 r = 1.5*_K/(_K+1)*cosf(2*M_PI*_fd*t*cosf(_theta));
448
449 // Window
450 w = kaiser(i, _n, beta, 0);
451
452 // composite
453 _h[i] = (J+r)*w;
454
455 //printf("t=%f, J=%f, r=%f, w=%f\n", t, J, r, w);
456 }
457 }
458
459
460 //
461 // filter analysis
462 //
463
464 // liquid_filter_autocorr()
465 //
466 // Compute auto-correlation of filter at a specific lag.
467 //
468 // _h : filter coefficients [size: _h_len x 1]
469 // _h_len : filter length
470 // _lag : auto-correlation lag (samples)
liquid_filter_autocorr(float * _h,unsigned int _h_len,int _lag)471 float liquid_filter_autocorr(float * _h,
472 unsigned int _h_len,
473 int _lag)
474 {
475 // auto-correlation is even symmetric
476 _lag = abs(_lag);
477
478 // lag outside of filter length is zero
479 if (_lag >= _h_len) return 0.0f;
480
481 // compute auto-correlation
482 float rxx=0.0f; // initialize auto-correlation to zero
483 unsigned int i;
484 for (i=_lag; i<_h_len; i++)
485 rxx += _h[i] * _h[i-_lag];
486
487 return rxx;
488 }
489
490 // liquid_filter_crosscorr()
491 //
492 // Compute cross-correlation of two filters at a specific lag.
493 //
494 // _h : filter coefficients [size: _h_len]
495 // _h_len : filter length
496 // _g : filter coefficients [size: _g_len]
497 // _g_len : filter length
498 // _lag : cross-correlation lag (samples)
liquid_filter_crosscorr(float * _h,unsigned int _h_len,float * _g,unsigned int _g_len,int _lag)499 float liquid_filter_crosscorr(float * _h,
500 unsigned int _h_len,
501 float * _g,
502 unsigned int _g_len,
503 int _lag)
504 {
505 // cross-correlation is odd symmetric
506 if (_h_len < _g_len) {
507 return liquid_filter_crosscorr(_g, _g_len,
508 _h, _h_len,
509 -_lag);
510 }
511
512 // at this point _h_len > _g_len
513 // assert(_h_len > _g_len);
514
515 if (_lag <= -(int)_g_len) return 0.0f;
516 if (_lag >= (int)_h_len) return 0.0f;
517
518 int ig = _lag < 0 ? -_lag : 0; // starting index for _g
519 int ih = _lag > 0 ? _lag : 0; // starting index for _h
520
521 // compute length of overlap
522 // condition 1: condition 2: condition 3:
523 // [------ h ------] [------ h ------] [------ h ------]
524 // [-- g --] [-- g --] [-- g --]
525 // >| n |< >| n |< >| n |<
526 //
527 int n;
528 if (_lag < 0)
529 n = (int)_g_len + _lag;
530 else if (_lag < (_h_len-_g_len))
531 n = _g_len;
532 else
533 n = _h_len - _lag;
534
535 // compute cross-correlation
536 float rxy=0.0f; // initialize auto-correlation to zero
537 int i;
538 for (i=0; i< n; i++)
539 rxy += _h[ih+i] * _g[ig+i];
540
541 return rxy;
542 }
543
544 // liquid_filter_isi()
545 //
546 // Compute inter-symbol interference (ISI)--both RMS and
547 // maximum--for the filter _h.
548 //
549 // _h : filter coefficients [size: 2*_k*_m+1 x 1]
550 // _k : filter over-sampling rate (samples/symbol)
551 // _m : filter delay (symbols)
552 // _rms : output root mean-squared ISI
553 // _max : maximum ISI
liquid_filter_isi(float * _h,unsigned int _k,unsigned int _m,float * _rms,float * _max)554 void liquid_filter_isi(float * _h,
555 unsigned int _k,
556 unsigned int _m,
557 float * _rms,
558 float * _max)
559 {
560 unsigned int h_len = 2*_k*_m+1;
561
562 // compute zero-lag auto-correlation
563 float rxx0 = liquid_filter_autocorr(_h,h_len,0);
564 //printf("rxx0 = %12.8f\n", rxx0);
565 //exit(1);
566
567 unsigned int i;
568 float isi_rms = 0.0f;
569 float isi_max = 0.0f;
570 float e;
571 for (i=1; i<=2*_m; i++) {
572 e = liquid_filter_autocorr(_h,h_len,i*_k) / rxx0;
573 e = fabsf(e);
574
575 isi_rms += e*e;
576
577 if (i==1 || e > isi_max)
578 isi_max = e;
579 }
580
581 *_rms = sqrtf( isi_rms / (float)(2*_m) );
582 *_max = isi_max;
583 }
584
585 // Compute relative out-of-band energy
586 //
587 // _h : filter coefficients [size: _h_len x 1]
588 // _h_len : filter length
589 // _fc : analysis cut-off frequency
590 // _nfft : fft size
liquid_filter_energy(float * _h,unsigned int _h_len,float _fc,unsigned int _nfft)591 float liquid_filter_energy(float * _h,
592 unsigned int _h_len,
593 float _fc,
594 unsigned int _nfft)
595 {
596 // validate input
597 if (_fc < 0.0f || _fc > 0.5f) {
598 fprintf(stderr,"error: liquid_filter_energy(), cut-off frequency must be in [0,0.5]\n");
599 exit(1);
600 } else if (_h_len == 0) {
601 fprintf(stderr,"error: liquid_filter_energy(), filter length must be greater than zero\n");
602 exit(1);
603 } else if (_nfft == 0) {
604 fprintf(stderr,"error: liquid_filter_energy(), fft size must be greater than zero\n");
605 exit(1);
606 }
607
608 // allocate memory for complex phasor
609 float complex expjwt[_h_len];
610
611 // initialize accumulators
612 float e_total = 0.0f; // total energy
613 float e_stopband = 0.0f; // stop-band energy
614
615 // create dotprod object
616 dotprod_crcf dp = dotprod_crcf_create(_h,_h_len);
617
618 unsigned int i;
619 unsigned int k;
620 for (i=0; i<_nfft; i++) {
621 float f = 0.5f * (float)i / (float)(_nfft);
622
623 // initialize complex phasor
624 for (k=0; k<_h_len; k++)
625 expjwt[k] = cexpf(_Complex_I*2*M_PI*f*k);
626
627 // compute vector dot product
628 float complex v;
629 dotprod_crcf_execute(dp, expjwt, &v);
630
631 // accumulate output
632 float e2 = crealf( v*conjf(v) );
633 e_total += e2;
634 e_stopband += (f >= _fc) ? e2 : 0.0f;
635 }
636
637 // destroy dotprod object
638 dotprod_crcf_destroy(dp);
639
640 // return energy ratio
641 return e_stopband / e_total;
642 }
643
644 // returns filter type based on input string
liquid_getopt_str2firfilt(const char * _str)645 int liquid_getopt_str2firfilt(const char * _str)
646 {
647 // Generic filter designs
648 if (strcmp(_str,"kaiser") ==0) return LIQUID_FIRFILT_KAISER;
649 else if (strcmp(_str,"pm") ==0) return LIQUID_FIRFILT_PM;
650
651 // Nyquist filter designs
652 else if (strcmp(_str,"rcos") ==0) return LIQUID_FIRFILT_RCOS;
653 else if (strcmp(_str,"fexp") ==0) return LIQUID_FIRFILT_FEXP;
654 else if (strcmp(_str,"fsech") ==0) return LIQUID_FIRFILT_FSECH;
655 else if (strcmp(_str,"farcsech") ==0) return LIQUID_FIRFILT_FARCSECH;
656
657 // root-Nyquist filter designs
658 else if (strcmp(_str,"arkaiser") ==0) return LIQUID_FIRFILT_ARKAISER;
659 else if (strcmp(_str,"rkaiser") ==0) return LIQUID_FIRFILT_RKAISER;
660 else if (strcmp(_str,"rrcos") ==0) return LIQUID_FIRFILT_RRC;
661 else if (strcmp(_str,"hM3") ==0) return LIQUID_FIRFILT_hM3;
662 else if (strcmp(_str,"gmsktx") ==0) return LIQUID_FIRFILT_GMSKTX;
663 else if (strcmp(_str,"gmskrx") ==0) return LIQUID_FIRFILT_GMSKRX;
664 else if (strcmp(_str,"rfexp") ==0) return LIQUID_FIRFILT_RFEXP;
665 else if (strcmp(_str,"rfsech") ==0) return LIQUID_FIRFILT_RFSECH;
666 else if (strcmp(_str,"rfarcsech")==0) return LIQUID_FIRFILT_RFARCSECH;
667
668 // filter type unknown
669 return LIQUID_FIRFILT_UNKNOWN;
670 }
671
672
673
674