3 * Copyright 2002,2007,2008 Free Software Foundation, Inc.
5 * This file is part of GNU Radio
7 * GNU Radio is free software; you can redistribute it and/or modify
8 * it under the terms of the GNU General Public License as published by
9 * the Free Software Foundation; either version 3, or (at your option)
12 * GNU Radio is distributed in the hope that it will be useful,
13 * but WITHOUT ANY WARRANTY; without even the implied warranty of
14 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
15 * GNU General Public License for more details.
17 * You should have received a copy of the GNU General Public License
18 * along with GNU Radio; see the file COPYING. If not, write to
19 * the Free Software Foundation, Inc., 51 Franklin Street,
20 * Boston, MA 02110-1301, USA.
23 #include <gr_firdes.h>
29 #define IzeroEPSILON 1E-21 /* Max error acceptable in Izero */
31 static double Izero(double x)
33 double sum, u, halfx, temp;
39 temp = halfx/(double)n;
44 } while (u >= IzeroEPSILON*sum);
54 gr_firdes::low_pass_2(double gain,
55 double sampling_freq, // Hz
56 double cutoff_freq, // Hz BEGINNING of transition band
57 double transition_width, // Hz width of transition band
58 double attenuation_dB, // attenuation dB
60 double beta) // used only with Kaiser
62 sanity_check_1f (sampling_freq, cutoff_freq, transition_width);
64 int ntaps = compute_ntaps_windes (sampling_freq, transition_width,
67 // construct the truncated ideal impulse response
68 // [sin(x)/x for the low pass case]
70 vector<float> taps(ntaps);
71 vector<float> w = window (window_type, ntaps, beta);
73 int M = (ntaps - 1) / 2;
74 double fwT0 = 2 * M_PI * cutoff_freq / sampling_freq;
75 for (int n = -M; n <= M; n++){
77 taps[n + M] = fwT0 / M_PI * w[n + M];
79 // a little algebra gets this into the more familiar sin(x)/x form
80 taps[n + M] = sin (n * fwT0) / (n * M_PI) * w[n + M];
84 // find the factor to normalize the gain, fmax.
85 // For low-pass, gain @ zero freq = 1.0
87 double fmax = taps[0 + M];
88 for (int n = 1; n <= M; n++)
89 fmax += 2 * taps[n + M];
91 gain /= fmax; // normalize
93 for (int i = 0; i < ntaps; i++)
101 gr_firdes::low_pass (double gain,
102 double sampling_freq,
103 double cutoff_freq, // Hz center of transition band
104 double transition_width, // Hz width of transition band
105 win_type window_type,
106 double beta) // used only with Kaiser
108 sanity_check_1f (sampling_freq, cutoff_freq, transition_width);
110 int ntaps = compute_ntaps (sampling_freq, transition_width,
113 // construct the truncated ideal impulse response
114 // [sin(x)/x for the low pass case]
116 vector<float> taps(ntaps);
117 vector<float> w = window (window_type, ntaps, beta);
119 int M = (ntaps - 1) / 2;
120 double fwT0 = 2 * M_PI * cutoff_freq / sampling_freq;
122 for (int n = -M; n <= M; n++){
124 taps[n + M] = fwT0 / M_PI * w[n + M];
126 // a little algebra gets this into the more familiar sin(x)/x form
127 taps[n + M] = sin (n * fwT0) / (n * M_PI) * w[n + M];
131 // find the factor to normalize the gain, fmax.
132 // For low-pass, gain @ zero freq = 1.0
134 double fmax = taps[0 + M];
135 for (int n = 1; n <= M; n++)
136 fmax += 2 * taps[n + M];
138 gain /= fmax; // normalize
140 for (int i = 0; i < ntaps; i++)
152 gr_firdes::high_pass_2 (double gain,
153 double sampling_freq,
154 double cutoff_freq, // Hz center of transition band
155 double transition_width, // Hz width of transition band
156 double attenuation_dB, // attenuation dB
157 win_type window_type,
158 double beta) // used only with Kaiser
160 sanity_check_1f (sampling_freq, cutoff_freq, transition_width);
162 int ntaps = compute_ntaps_windes (sampling_freq, transition_width,
165 // construct the truncated ideal impulse response times the window function
167 vector<float> taps(ntaps);
168 vector<float> w = window (window_type, ntaps, beta);
170 int M = (ntaps - 1) / 2;
171 double fwT0 = 2 * M_PI * cutoff_freq / sampling_freq;
173 for (int n = -M; n <= M; n++){
175 taps[n + M] = (1 - (fwT0 / M_PI)) * w[n + M];
177 // a little algebra gets this into the more familiar sin(x)/x form
178 taps[n + M] = -sin (n * fwT0) / (n * M_PI) * w[n + M];
182 // find the factor to normalize the gain, fmax.
183 // For high-pass, gain @ fs/2 freq = 1.0
185 double fmax = taps[0 + M];
186 for (int n = 1; n <= M; n++)
187 fmax += 2 * taps[n + M] * cos (n * M_PI);
189 gain /= fmax; // normalize
191 for (int i = 0; i < ntaps; i++)
200 gr_firdes::high_pass (double gain,
201 double sampling_freq,
202 double cutoff_freq, // Hz center of transition band
203 double transition_width, // Hz width of transition band
204 win_type window_type,
205 double beta) // used only with Kaiser
207 sanity_check_1f (sampling_freq, cutoff_freq, transition_width);
209 int ntaps = compute_ntaps (sampling_freq, transition_width,
212 // construct the truncated ideal impulse response times the window function
214 vector<float> taps(ntaps);
215 vector<float> w = window (window_type, ntaps, beta);
217 int M = (ntaps - 1) / 2;
218 double fwT0 = 2 * M_PI * cutoff_freq / sampling_freq;
220 for (int n = -M; n <= M; n++){
222 taps[n + M] = (1 - (fwT0 / M_PI)) * w[n + M];
224 // a little algebra gets this into the more familiar sin(x)/x form
225 taps[n + M] = -sin (n * fwT0) / (n * M_PI) * w[n + M];
229 // find the factor to normalize the gain, fmax.
230 // For high-pass, gain @ fs/2 freq = 1.0
232 double fmax = taps[0 + M];
233 for (int n = 1; n <= M; n++)
234 fmax += 2 * taps[n + M] * cos (n * M_PI);
236 gain /= fmax; // normalize
238 for (int i = 0; i < ntaps; i++)
249 gr_firdes::band_pass_2 (double gain,
250 double sampling_freq,
251 double low_cutoff_freq, // Hz center of transition band
252 double high_cutoff_freq, // Hz center of transition band
253 double transition_width, // Hz width of transition band
254 double attenuation_dB, // attenuation dB
255 win_type window_type,
256 double beta) // used only with Kaiser
258 sanity_check_2f (sampling_freq,
260 high_cutoff_freq, transition_width);
262 int ntaps = compute_ntaps_windes (sampling_freq, transition_width,
265 vector<float> taps(ntaps);
266 vector<float> w = window (window_type, ntaps, beta);
268 int M = (ntaps - 1) / 2;
269 double fwT0 = 2 * M_PI * low_cutoff_freq / sampling_freq;
270 double fwT1 = 2 * M_PI * high_cutoff_freq / sampling_freq;
272 for (int n = -M; n <= M; n++){
274 taps[n + M] = (fwT1 - fwT0) / M_PI * w[n + M];
276 taps[n + M] = (sin (n * fwT1) - sin (n * fwT0)) / (n * M_PI) * w[n + M];
280 // find the factor to normalize the gain, fmax.
281 // For band-pass, gain @ center freq = 1.0
283 double fmax = taps[0 + M];
284 for (int n = 1; n <= M; n++)
285 fmax += 2 * taps[n + M] * cos (n * (fwT0 + fwT1) * 0.5);
287 gain /= fmax; // normalize
289 for (int i = 0; i < ntaps; i++)
297 gr_firdes::band_pass (double gain,
298 double sampling_freq,
299 double low_cutoff_freq, // Hz center of transition band
300 double high_cutoff_freq, // Hz center of transition band
301 double transition_width, // Hz width of transition band
302 win_type window_type,
303 double beta) // used only with Kaiser
305 sanity_check_2f (sampling_freq,
307 high_cutoff_freq, transition_width);
309 int ntaps = compute_ntaps (sampling_freq, transition_width,
312 // construct the truncated ideal impulse response times the window function
314 vector<float> taps(ntaps);
315 vector<float> w = window (window_type, ntaps, beta);
317 int M = (ntaps - 1) / 2;
318 double fwT0 = 2 * M_PI * low_cutoff_freq / sampling_freq;
319 double fwT1 = 2 * M_PI * high_cutoff_freq / sampling_freq;
321 for (int n = -M; n <= M; n++){
323 taps[n + M] = (fwT1 - fwT0) / M_PI * w[n + M];
325 taps[n + M] = (sin (n * fwT1) - sin (n * fwT0)) / (n * M_PI) * w[n + M];
329 // find the factor to normalize the gain, fmax.
330 // For band-pass, gain @ center freq = 1.0
332 double fmax = taps[0 + M];
333 for (int n = 1; n <= M; n++)
334 fmax += 2 * taps[n + M] * cos (n * (fwT0 + fwT1) * 0.5);
336 gain /= fmax; // normalize
338 for (int i = 0; i < ntaps; i++)
345 // === Complex Band Pass ===
349 gr_firdes::complex_band_pass_2 (double gain,
350 double sampling_freq,
351 double low_cutoff_freq, // Hz center of transition band
352 double high_cutoff_freq, // Hz center of transition band
353 double transition_width, // Hz width of transition band
354 double attenuation_dB, // attenuation dB
355 win_type window_type,
356 double beta) // used only with Kaiser
358 sanity_check_2f_c (sampling_freq,
360 high_cutoff_freq, transition_width);
362 int ntaps = compute_ntaps_windes (sampling_freq, transition_width,
367 vector<gr_complex> taps(ntaps);
368 vector<float> lptaps(ntaps);
369 vector<float> w = window (window_type, ntaps, beta);
371 lptaps = low_pass_2(gain,sampling_freq,(high_cutoff_freq - low_cutoff_freq)/2,transition_width,attenuation_dB,window_type,beta);
373 gr_complex *optr = &taps[0];
374 float *iptr = &lptaps[0];
375 float freq = M_PI * (high_cutoff_freq + low_cutoff_freq)/sampling_freq;
377 if (lptaps.size() & 01) {
378 phase = - freq * ( lptaps.size() >> 1 );
379 } else phase = - freq/2.0 * ((1 + 2*lptaps.size()) >> 1);
380 for(unsigned int i=0;i<lptaps.size();i++) {
381 *optr++ = gr_complex(*iptr * cos(phase),*iptr * sin(phase));
382 iptr++, phase += freq;
390 gr_firdes::complex_band_pass (double gain,
391 double sampling_freq,
392 double low_cutoff_freq, // Hz center of transition band
393 double high_cutoff_freq, // Hz center of transition band
394 double transition_width, // Hz width of transition band
395 win_type window_type,
396 double beta) // used only with Kaiser
398 sanity_check_2f_c (sampling_freq,
400 high_cutoff_freq, transition_width);
402 int ntaps = compute_ntaps (sampling_freq, transition_width,
405 // construct the truncated ideal impulse response times the window function
407 vector<gr_complex> taps(ntaps);
408 vector<float> lptaps(ntaps);
409 vector<float> w = window (window_type, ntaps, beta);
411 lptaps = low_pass(gain,sampling_freq,(high_cutoff_freq - low_cutoff_freq)/2,transition_width,window_type,beta);
413 gr_complex *optr = &taps[0];
414 float *iptr = &lptaps[0];
415 float freq = M_PI * (high_cutoff_freq + low_cutoff_freq)/sampling_freq;
417 if (lptaps.size() & 01) {
418 phase = - freq * ( lptaps.size() >> 1 );
419 } else phase = - freq/2.0 * ((1 + 2*lptaps.size()) >> 1);
420 for(unsigned int i=0;i<lptaps.size();i++) {
421 *optr++ = gr_complex(*iptr * cos(phase),*iptr * sin(phase));
422 iptr++, phase += freq;
429 // === Band Reject ===
433 gr_firdes::band_reject_2 (double gain,
434 double sampling_freq,
435 double low_cutoff_freq, // Hz center of transition band
436 double high_cutoff_freq, // Hz center of transition band
437 double transition_width, // Hz width of transition band
438 double attenuation_dB, // attenuation dB
439 win_type window_type,
440 double beta) // used only with Kaiser
442 sanity_check_2f (sampling_freq,
444 high_cutoff_freq, transition_width);
446 int ntaps = compute_ntaps_windes (sampling_freq, transition_width,
449 // construct the truncated ideal impulse response times the window function
451 vector<float> taps(ntaps);
452 vector<float> w = window (window_type, ntaps, beta);
454 int M = (ntaps - 1) / 2;
455 double fwT0 = 2 * M_PI * low_cutoff_freq / sampling_freq;
456 double fwT1 = 2 * M_PI * high_cutoff_freq / sampling_freq;
458 for (int n = -M; n <= M; n++){
460 taps[n + M] = 1.0 + ((fwT0 - fwT1) / M_PI * w[n + M]);
462 taps[n + M] = (sin (n * fwT0) - sin (n * fwT1)) / (n * M_PI) * w[n + M];
466 // find the factor to normalize the gain, fmax.
467 // For band-reject, gain @ zero freq = 1.0
469 double fmax = taps[0 + M];
470 for (int n = 1; n <= M; n++)
471 fmax += 2 * taps[n + M];
473 gain /= fmax; // normalize
475 for (int i = 0; i < ntaps; i++)
482 gr_firdes::band_reject (double gain,
483 double sampling_freq,
484 double low_cutoff_freq, // Hz center of transition band
485 double high_cutoff_freq, // Hz center of transition band
486 double transition_width, // Hz width of transition band
487 win_type window_type,
488 double beta) // used only with Kaiser
490 sanity_check_2f (sampling_freq,
492 high_cutoff_freq, transition_width);
494 int ntaps = compute_ntaps (sampling_freq, transition_width,
497 // construct the truncated ideal impulse response times the window function
499 vector<float> taps(ntaps);
500 vector<float> w = window (window_type, ntaps, beta);
502 int M = (ntaps - 1) / 2;
503 double fwT0 = 2 * M_PI * low_cutoff_freq / sampling_freq;
504 double fwT1 = 2 * M_PI * high_cutoff_freq / sampling_freq;
506 for (int n = -M; n <= M; n++){
508 taps[n + M] = 1.0 + ((fwT0 - fwT1) / M_PI * w[n + M]);
510 taps[n + M] = (sin (n * fwT0) - sin (n * fwT1)) / (n * M_PI) * w[n + M];
514 // find the factor to normalize the gain, fmax.
515 // For band-reject, gain @ zero freq = 1.0
517 double fmax = taps[0 + M];
518 for (int n = 1; n <= M; n++)
519 fmax += 2 * taps[n + M];
521 gain /= fmax; // normalize
523 for (int i = 0; i < ntaps; i++)
534 gr_firdes::hilbert (unsigned int ntaps,
539 throw std::out_of_range("Hilbert: Must have odd number of taps");
541 vector<float> taps(ntaps);
542 vector<float> w = window (windowtype, ntaps, beta);
543 unsigned int h = (ntaps-1)/2;
545 for (unsigned int i = 1; i <= h; i++)
549 float x = 1/(float)i;
550 taps[h+i] = x * w[h+i];
551 taps[h-i] = -x * w[h-i];
552 gain = taps[h+i] - gain;
555 taps[h+i] = taps[h-i] = 0;
557 gain = 2 * fabs(gain);
558 for ( unsigned int i = 0; i < ntaps; i++)
568 gr_firdes::gaussian (double gain,
574 vector<float> taps(ntaps);
577 double s = 1.0/(sqrt(log(2)) / (2*M_PI*bt));
578 double t0 = -0.5 * ntaps;
580 for(int i=0;i<ntaps;i++)
584 taps[i] = exp(-0.5*ts*ts);
587 for(int i=0;i<ntaps;i++)
588 taps[i] = taps[i] / scale * gain;
594 // Root Raised Cosine
598 gr_firdes::root_raised_cosine (double gain,
599 double sampling_freq,
604 ntaps |= 1; // ensure that ntaps is odd
606 double spb = sampling_freq/symbol_rate; // samples per bit/symbol
607 vector<float> taps(ntaps);
609 for(int i=0;i<ntaps;i++)
611 double x1,x2,x3,num,den;
612 double xindx = i - ntaps/2;
613 x1 = M_PI * xindx/spb;
614 x2 = 4 * alpha * xindx / spb;
616 if( fabs(x3) >= 0.000001 ) // Avoid Rounding errors...
619 num = cos((1+alpha)*x1) + sin((1-alpha)*x1)/(4*alpha*xindx/spb);
621 num = cos((1+alpha)*x1) + (1-alpha) * M_PI / (4*alpha);
633 num = (sin(x2)*(1+alpha)*M_PI
634 - cos(x3)*((1-alpha)*M_PI*spb)/(4*alpha*xindx)
635 + sin(x3)*spb*spb/(4*alpha*xindx*xindx));
636 den = -32 * M_PI * alpha * alpha * xindx/spb;
638 taps[i] = 4 * alpha * num / den;
642 for(int i=0;i<ntaps;i++)
643 taps[i] = taps[i] * gain / scale;
652 // delta_f / width_factor gives number of taps required.
653 static const float width_factor[5] = { // indexed by win_type
657 2.0, // WIN_RECTANGULAR
658 //5.0 // WIN_KAISER (guesstimate compromise)
659 //2.0 // WIN_KAISER (guesstimate compromise)
664 gr_firdes::compute_ntaps_windes(double sampling_freq,
665 double transition_width, // this is frequency, not relative frequency
666 double attenuation_dB)
668 // Based on formula from Multirate Signal Processing for
669 // Communications Systems, fredric j harris
670 int ntaps = (int)(attenuation_dB*sampling_freq/(22.0*transition_width));
671 if ((ntaps & 1) == 0) // if even...
672 ntaps++; // ...make odd
677 gr_firdes::compute_ntaps (double sampling_freq,
678 double transition_width,
679 win_type window_type,
682 // normalized transition width
683 double delta_f = transition_width / sampling_freq;
685 // compute number of taps required for given transition width
686 int ntaps = (int) (width_factor[window_type] / delta_f + 0.5);
687 if ((ntaps & 1) == 0) // if even...
688 ntaps++; // ...make odd
693 double gr_firdes::bessi0(double x)
703 ans=1.0+y*(3.5156229+y*(3.0899424+y*(1.2067492
704 +y*(0.2659732+y*(0.360768e-1+y*0.45813e-2)))));
709 ans=(exp(ax)/sqrt(ax))*(0.39894228+y*(0.1328592e-1
710 +y*(0.225319e-2+y*(-0.157565e-2+y*(0.916281e-2
711 +y*(-0.2057706e-1+y*(0.2635537e-1+y*(-0.1647633e-1
712 +y*0.392377e-2))))))));
717 gr_firdes::window (win_type type, int ntaps, double beta)
719 vector<float> taps(ntaps);
720 int M = ntaps - 1; // filter order
723 case WIN_RECTANGULAR:
724 for (int n = 0; n < ntaps; n++)
728 for (int n = 0; n < ntaps; n++)
729 taps[n] = 0.54 - 0.46 * cos ((2 * M_PI * n) / M);
733 for (int n = 0; n < ntaps; n++)
734 taps[n] = 0.5 - 0.5 * cos ((2 * M_PI * n) / M);
738 for (int n = 0; n < ntaps; n++)
739 taps[n] = 0.42 - 0.50 * cos ((2*M_PI * n) / (M-1)) - 0.08 * cos ((4*M_PI * n) / (M-1));
742 case WIN_BLACKMAN_hARRIS:
743 for (int n = -ntaps/2; n < ntaps/2; n++)
744 taps[n+ntaps/2] = 0.35875 + 0.48829*cos((2*M_PI * n) / (float)M) +
745 0.14128*cos((4*M_PI * n) / (float)M) + 0.01168*cos((6*M_PI * n) / (float)M);
750 for (int n = 0; n < ntaps; n++)
751 taps[n] = bessi0(beta*sqrt(1.0 - (4.0*n/(M*M))))/bessi0(beta);
757 double IBeta = 1.0/Izero(beta);
758 double inm1 = 1.0/((double)(ntaps));
760 //fprintf(stderr, "IBeta = %g; inm1 = %g\n", IBeta, inm1);
762 for (int i=0; i<ntaps; i++) {
764 //fprintf(stderr, "temp = %g\n", temp);
765 taps[i] = Izero(beta*sqrt(1.0-temp*temp)) * IBeta;
766 //fprintf(stderr, "taps[%d] = %g\n", i, taps[i]);
773 throw std::out_of_range ("gr_firdes:window: type out of range");
780 gr_firdes::sanity_check_1f (double sampling_freq,
781 double fa, // cutoff freq
782 double transition_width)
784 if (sampling_freq <= 0.0)
785 throw std::out_of_range ("gr_firdes check failed: sampling_freq > 0");
787 if (fa <= 0.0 || fa > sampling_freq / 2)
788 throw std::out_of_range ("gr_firdes check failed: 0 < fa <= sampling_freq / 2");
790 if (transition_width <= 0)
791 throw std::out_of_range ("gr_dirdes check failed: transition_width > 0");
795 gr_firdes::sanity_check_2f (double sampling_freq,
796 double fa, // first cutoff freq
797 double fb, // second cutoff freq
798 double transition_width)
800 if (sampling_freq <= 0.0)
801 throw std::out_of_range ("gr_firdes check failed: sampling_freq > 0");
803 if (fa <= 0.0 || fa > sampling_freq / 2)
804 throw std::out_of_range ("gr_firdes check failed: 0 < fa <= sampling_freq / 2");
806 if (fb <= 0.0 || fb > sampling_freq / 2)
807 throw std::out_of_range ("gr_firdes check failed: 0 < fb <= sampling_freq / 2");
810 throw std::out_of_range ("gr_firdes check failed: fa <= fb");
812 if (transition_width <= 0)
813 throw std::out_of_range ("gr_firdes check failed: transition_width > 0");
817 gr_firdes::sanity_check_2f_c (double sampling_freq,
818 double fa, // first cutoff freq
819 double fb, // second cutoff freq
820 double transition_width)
822 if (sampling_freq <= 0.0)
823 throw std::out_of_range ("gr_firdes check failed: sampling_freq > 0");
825 if (fa < -sampling_freq / 2 || fa > sampling_freq / 2)
826 throw std::out_of_range ("gr_firdes check failed: 0 < fa <= sampling_freq / 2");
828 if (fb < -sampling_freq / 2 || fb > sampling_freq / 2)
829 throw std::out_of_range ("gr_firdes check failed: 0 < fb <= sampling_freq / 2");
832 throw std::out_of_range ("gr_firdes check failed: fa <= fb");
834 if (transition_width <= 0)
835 throw std::out_of_range ("gr_firdes check failed: transition_width > 0");