3 * Copyright 2002 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 2, 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 (double gain,
56 double cutoff_freq, // Hz center of transition band
57 double transition_width, // Hz width of transition band
59 double beta) // used only with Kaiser
61 sanity_check_1f (sampling_freq, cutoff_freq, transition_width);
63 int ntaps = compute_ntaps (sampling_freq, transition_width,
66 // construct the truncated ideal impulse response
67 // [sin(x)/x for the low pass case]
69 vector<float> taps(ntaps);
70 vector<float> w = window (window_type, ntaps, beta);
72 int M = (ntaps - 1) / 2;
73 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++)
104 gr_firdes::high_pass (double gain,
105 double sampling_freq,
106 double cutoff_freq, // Hz center of transition band
107 double transition_width, // Hz width of transition band
108 win_type window_type,
109 double beta) // used only with Kaiser
111 sanity_check_1f (sampling_freq, cutoff_freq, transition_width);
113 int ntaps = compute_ntaps (sampling_freq, transition_width,
116 // construct the truncated ideal impulse response times the window function
118 vector<float> taps(ntaps);
119 vector<float> w = window (window_type, ntaps, beta);
121 int M = (ntaps - 1) / 2;
122 double fwT0 = 2 * M_PI * cutoff_freq / sampling_freq;
124 for (int n = -M; n <= M; n++){
126 taps[n + M] = (1 - (fwT0 / M_PI)) * w[n + M];
128 // a little algebra gets this into the more familiar sin(x)/x form
129 taps[n + M] = -sin (n * fwT0) / (n * M_PI) * w[n + M];
133 // find the factor to normalize the gain, fmax.
134 // For high-pass, gain @ fs/2 freq = 1.0
136 double fmax = taps[0 + M];
137 for (int n = 1; n <= M; n++)
138 fmax += 2 * taps[n + M] * cos (n * M_PI);
140 gain /= fmax; // normalize
142 for (int i = 0; i < ntaps; i++)
153 gr_firdes::band_pass (double gain,
154 double sampling_freq,
155 double low_cutoff_freq, // Hz center of transition band
156 double high_cutoff_freq, // Hz center of transition band
157 double transition_width, // Hz width of transition band
158 win_type window_type,
159 double beta) // used only with Kaiser
161 sanity_check_2f (sampling_freq,
163 high_cutoff_freq, transition_width);
165 int ntaps = compute_ntaps (sampling_freq, transition_width,
168 // construct the truncated ideal impulse response times the window function
170 vector<float> taps(ntaps);
171 vector<float> w = window (window_type, ntaps, beta);
173 int M = (ntaps - 1) / 2;
174 double fwT0 = 2 * M_PI * low_cutoff_freq / sampling_freq;
175 double fwT1 = 2 * M_PI * high_cutoff_freq / sampling_freq;
177 for (int n = -M; n <= M; n++){
179 taps[n + M] = (fwT1 - fwT0) / M_PI * w[n + M];
181 taps[n + M] = (sin (n * fwT1) - sin (n * fwT0)) / (n * M_PI) * w[n + M];
185 // find the factor to normalize the gain, fmax.
186 // For band-pass, gain @ center freq = 1.0
188 double fmax = taps[0 + M];
189 for (int n = 1; n <= M; n++)
190 fmax += 2 * taps[n + M] * cos (n * (fwT0 + fwT1) * 0.5);
192 gain /= fmax; // normalize
194 for (int i = 0; i < ntaps; i++)
201 // === Complex Band Pass ===
205 gr_firdes::complex_band_pass (double gain,
206 double sampling_freq,
207 double low_cutoff_freq, // Hz center of transition band
208 double high_cutoff_freq, // Hz center of transition band
209 double transition_width, // Hz width of transition band
210 win_type window_type,
211 double beta) // used only with Kaiser
213 sanity_check_2f_c (sampling_freq,
215 high_cutoff_freq, transition_width);
217 int ntaps = compute_ntaps (sampling_freq, transition_width,
220 // construct the truncated ideal impulse response times the window function
222 vector<gr_complex> taps(ntaps);
223 vector<float> lptaps(ntaps);
224 vector<float> w = window (window_type, ntaps, beta);
226 lptaps = low_pass(gain,sampling_freq,(high_cutoff_freq - low_cutoff_freq)/2,transition_width,window_type,beta);
228 gr_complex *optr = &taps[0];
229 float *iptr = &lptaps[0];
230 float freq = M_PI * (high_cutoff_freq + low_cutoff_freq)/sampling_freq;
232 if (lptaps.size() & 01) {
233 phase = - freq * ( lptaps.size() >> 1 );
234 } else phase = - freq/2.0 * ((1 + 2*lptaps.size()) >> 1);
235 for(unsigned int i=0;i<lptaps.size();i++) {
236 *optr++ = gr_complex(*iptr * cos(phase),*iptr * sin(phase));
237 iptr++, phase += freq;
245 // === Band Reject ===
249 gr_firdes::band_reject (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 win_type window_type,
255 double beta) // used only with Kaiser
257 sanity_check_2f (sampling_freq,
259 high_cutoff_freq, transition_width);
261 int ntaps = compute_ntaps (sampling_freq, transition_width,
264 // construct the truncated ideal impulse response times the window function
266 vector<float> taps(ntaps);
267 vector<float> w = window (window_type, ntaps, beta);
269 int M = (ntaps - 1) / 2;
270 double fwT0 = 2 * M_PI * low_cutoff_freq / sampling_freq;
271 double fwT1 = 2 * M_PI * high_cutoff_freq / sampling_freq;
273 for (int n = -M; n <= M; n++){
275 taps[n + M] = (1.0 + (fwT0 - fwT1)) / M_PI * w[n + M];
277 taps[n + M] = (sin (n * fwT0) - sin (n * fwT1)) / (n * M_PI) * w[n + M];
281 // find the factor to normalize the gain, fmax.
282 // For band-reject, gain @ zero freq = 1.0
284 double fmax = taps[0 + M];
285 for (int n = 1; n <= M; n++)
286 fmax += 2 * taps[n + M];
288 gain /= fmax; // normalize
290 for (int i = 0; i < ntaps; i++)
301 gr_firdes::hilbert (unsigned int ntaps,
306 throw std::out_of_range("Hilbert: Must have odd number of taps");
308 vector<float> taps(ntaps);
309 vector<float> w = window (windowtype, ntaps, beta);
310 unsigned int h = (ntaps-1)/2;
312 for (unsigned int i = 1; i <= h; i++)
316 float x = 1/(float)i;
317 taps[h+i] = x * w[h+i];
318 taps[h-i] = -x * w[h-i];
319 gain = taps[h+i] - gain;
322 taps[h+i] = taps[h-i] = 0;
324 gain = 2 * fabs(gain);
325 for ( unsigned int i = 0; i < ntaps; i++)
335 gr_firdes::gaussian (double gain,
341 vector<float> taps(ntaps);
344 double s = 1.0/(sqrt(log(2)) / (2*M_PI*bt));
345 double t0 = -0.5 * ntaps;
347 for(int i=0;i<ntaps;i++)
351 taps[i] = exp(-0.5*ts*ts);
354 for(int i=0;i<ntaps;i++)
355 taps[i] = taps[i] / scale * gain;
361 // Root Raised Cosine
365 gr_firdes::root_raised_cosine (double gain,
366 double sampling_freq,
371 ntaps |= 1; // ensure that ntaps is odd
373 double spb = sampling_freq/symbol_rate; // samples per bit/symbol
374 vector<float> taps(ntaps);
376 for(int i=0;i<ntaps;i++)
378 double x1,x2,x3,num,den;
379 double xindx = i - ntaps/2;
380 x1 = M_PI * xindx/spb;
381 x2 = 4 * alpha * xindx / spb;
383 if( fabs(x3) >= 0.000001 ) // Avoid Rounding errors...
386 num = cos((1+alpha)*x1) + sin((1-alpha)*x1)/(4*alpha*xindx/spb);
388 num = cos((1+alpha)*x1) + (1-alpha) * M_PI / (4*alpha);
400 num = (sin(x2)*(1+alpha)*M_PI
401 - cos(x3)*((1-alpha)*M_PI*spb)/(4*alpha*xindx)
402 + sin(x3)*spb*spb/(4*alpha*xindx*xindx));
403 den = -32 * M_PI * alpha * alpha * xindx/spb;
405 taps[i] = 4 * alpha * num / den;
409 for(int i=0;i<ntaps;i++)
410 taps[i] = taps[i] * gain / scale;
419 // delta_f / width_factor gives number of taps required.
420 static const float width_factor[5] = { // indexed by win_type
424 2.0, // WIN_RECTANGULAR
425 //5.0 // WIN_KAISER (guesstimate compromise)
426 //2.0 // WIN_KAISER (guesstimate compromise)
431 gr_firdes::compute_ntaps (double sampling_freq,
432 double transition_width,
433 win_type window_type,
436 // normalized transition width
437 double delta_f = transition_width / sampling_freq;
439 // compute number of taps required for given transition width
440 int ntaps = (int) (width_factor[window_type] / delta_f + 0.5);
441 if ((ntaps & 1) == 0) // if even...
442 ntaps++; // ...make odd
447 double gr_firdes::bessi0(double x)
457 ans=1.0+y*(3.5156229+y*(3.0899424+y*(1.2067492
458 +y*(0.2659732+y*(0.360768e-1+y*0.45813e-2)))));
463 ans=(exp(ax)/sqrt(ax))*(0.39894228+y*(0.1328592e-1
464 +y*(0.225319e-2+y*(-0.157565e-2+y*(0.916281e-2
465 +y*(-0.2057706e-1+y*(0.2635537e-1+y*(-0.1647633e-1
466 +y*0.392377e-2))))))));
471 gr_firdes::window (win_type type, int ntaps, double beta)
473 vector<float> taps(ntaps);
474 int M = ntaps - 1; // filter order
477 case WIN_RECTANGULAR:
478 for (int n = 0; n < ntaps; n++)
482 for (int n = 0; n < ntaps; n++)
483 taps[n] = 0.54 - 0.46 * cos ((2 * M_PI * n) / M);
487 for (int n = 0; n < ntaps; n++)
488 taps[n] = 0.5 - 0.5 * cos ((2 * M_PI * n) / M);
492 for (int n = 0; n < ntaps; n++)
493 taps[n] = 0.42 - 0.50 * cos ((2*M_PI * n) / (M-1)) - 0.08 * cos ((4*M_PI * n) / (M-1));
498 for (int n = 0; n < ntaps; n++)
499 taps[n] = bessi0(beta*sqrt(1.0 - (4.0*n/(M*M))))/bessi0(beta);
505 double IBeta = 1.0/Izero(beta);
506 double inm1 = 1.0/((double)(ntaps));
508 //fprintf(stderr, "IBeta = %g; inm1 = %g\n", IBeta, inm1);
510 for (int i=0; i<ntaps; i++) {
512 //fprintf(stderr, "temp = %g\n", temp);
513 taps[i] = Izero(beta*sqrt(1.0-temp*temp)) * IBeta;
514 //fprintf(stderr, "taps[%d] = %g\n", i, taps[i]);
521 throw std::runtime_error ("not_implemented");
528 gr_firdes::sanity_check_1f (double sampling_freq,
529 double fa, // cutoff freq
530 double transition_width)
532 if (sampling_freq <= 0.0)
533 throw std::out_of_range ("gr_firdes check failed: sampling_freq > 0");
535 if (fa <= 0.0 || fa > sampling_freq / 2)
536 throw std::out_of_range ("gr_firdes check failed: 0 < fa <= sampling_freq / 2");
538 if (transition_width <= 0)
539 throw std::out_of_range ("gr_dirdes check failed: transition_width > 0");
543 gr_firdes::sanity_check_2f (double sampling_freq,
544 double fa, // first cutoff freq
545 double fb, // second cutoff freq
546 double transition_width)
548 if (sampling_freq <= 0.0)
549 throw std::out_of_range ("gr_firdes check failed: sampling_freq > 0");
551 if (fa <= 0.0 || fa > sampling_freq / 2)
552 throw std::out_of_range ("gr_firdes check failed: 0 < fa <= sampling_freq / 2");
554 if (fb <= 0.0 || fb > sampling_freq / 2)
555 throw std::out_of_range ("gr_firdes check failed: 0 < fb <= sampling_freq / 2");
558 throw std::out_of_range ("gr_firdes check failed: fa <= fb");
560 if (transition_width <= 0)
561 throw std::out_of_range ("gr_firdes check failed: transition_width > 0");
565 gr_firdes::sanity_check_2f_c (double sampling_freq,
566 double fa, // first cutoff freq
567 double fb, // second cutoff freq
568 double transition_width)
570 if (sampling_freq <= 0.0)
571 throw std::out_of_range ("gr_firdes check failed: sampling_freq > 0");
573 if (fa < -sampling_freq / 2 || fa > sampling_freq / 2)
574 throw std::out_of_range ("gr_firdes check failed: 0 < fa <= sampling_freq / 2");
576 if (fb < -sampling_freq / 2 || fb > sampling_freq / 2)
577 throw std::out_of_range ("gr_firdes check failed: 0 < fb <= sampling_freq / 2");
580 throw std::out_of_range ("gr_firdes check failed: fa <= fb");
582 if (transition_width <= 0)
583 throw std::out_of_range ("gr_firdes check failed: transition_width > 0");