From a3c5490e99dc21f24432c984cbdf0763e585a391 Mon Sep 17 00:00:00 2001 From: n4hy Date: Mon, 18 Aug 2008 22:36:19 +0000 Subject: New taps computation based on requested transition bandwidth and stopband attenuation. qa code testing each added git-svn-id: http://gnuradio.org/svn/gnuradio/trunk@9316 221aa14e-8319-0410-a670-987f0aec2ac5 --- gnuradio-core/src/lib/general/gr_firdes.cc | 270 +++++++++++++++++++++++-- gnuradio-core/src/lib/general/gr_firdes.h | 51 +++++ gnuradio-core/src/lib/general/qa_gr_firdes.cc | 272 ++++++++++++++++++++++++++ gnuradio-core/src/lib/general/qa_gr_firdes.h | 6 + 4 files changed, 587 insertions(+), 12 deletions(-) (limited to 'gnuradio-core/src/lib/general') diff --git a/gnuradio-core/src/lib/general/gr_firdes.cc b/gnuradio-core/src/lib/general/gr_firdes.cc index 06730b7f8..17c40d78a 100644 --- a/gnuradio-core/src/lib/general/gr_firdes.cc +++ b/gnuradio-core/src/lib/general/gr_firdes.cc @@ -50,6 +50,53 @@ static double Izero(double x) // === Low Pass === // +vector +gr_firdes::low_pass_2(double gain, + double sampling_freq, // Hz + double cutoff_freq, // Hz BEGINNING of transition band + double transition_width, // Hz width of transition band + double attenuation_dB, // attenuation dB + win_type window_type, + double beta) // used only with Kaiser +{ + sanity_check_1f (sampling_freq, cutoff_freq, transition_width); + + int ntaps = compute_ntaps_windes (sampling_freq, transition_width, + attenuation_dB); + + // construct the truncated ideal impulse response + // [sin(x)/x for the low pass case] + + vector taps(ntaps); + vector w = window (window_type, ntaps, beta); + + int M = (ntaps - 1) / 2; + double fwT0 = 2 * M_PI * cutoff_freq / sampling_freq; + for (int n = -M; n <= M; n++){ + if (n == 0) + taps[n + M] = fwT0 / M_PI * w[n + M]; + else { + // a little algebra gets this into the more familiar sin(x)/x form + taps[n + M] = sin (n * fwT0) / (n * M_PI) * w[n + M]; + } + } + + // find the factor to normalize the gain, fmax. + // For low-pass, gain @ zero freq = 1.0 + + double fmax = taps[0 + M]; + for (int n = 1; n <= M; n++) + fmax += 2 * taps[n + M]; + + gain /= fmax; // normalize + + for (int i = 0; i < ntaps; i++) + taps[i] *= gain; + + + return taps; +} + vector gr_firdes::low_pass (double gain, double sampling_freq, @@ -74,7 +121,7 @@ gr_firdes::low_pass (double gain, for (int n = -M; n <= M; n++){ if (n == 0) - taps[n + M] = fwT0 / M_PI * w[n + M]; + taps[n + M] = fwT0 / M_PI * w[n + M]; else { // a little algebra gets this into the more familiar sin(x)/x form taps[n + M] = sin (n * fwT0) / (n * M_PI) * w[n + M]; @@ -96,22 +143,71 @@ gr_firdes::low_pass (double gain, return taps; } + // // === High Pass === // +vector +gr_firdes::high_pass_2 (double gain, + double sampling_freq, + double cutoff_freq, // Hz center of transition band + double transition_width, // Hz width of transition band + double attenuation_dB, // attenuation dB + win_type window_type, + double beta) // used only with Kaiser +{ + sanity_check_1f (sampling_freq, cutoff_freq, transition_width); + + int ntaps = compute_ntaps_windes (sampling_freq, transition_width, + attenuation_dB); + + // construct the truncated ideal impulse response times the window function + + vector taps(ntaps); + vector w = window (window_type, ntaps, beta); + + int M = (ntaps - 1) / 2; + double fwT0 = 2 * M_PI * cutoff_freq / sampling_freq; + + for (int n = -M; n <= M; n++){ + if (n == 0) + taps[n + M] = (1 - (fwT0 / M_PI)) * w[n + M]; + else { + // a little algebra gets this into the more familiar sin(x)/x form + taps[n + M] = -sin (n * fwT0) / (n * M_PI) * w[n + M]; + } + } + + // find the factor to normalize the gain, fmax. + // For high-pass, gain @ fs/2 freq = 1.0 + + double fmax = taps[0 + M]; + for (int n = 1; n <= M; n++) + fmax += 2 * taps[n + M] * cos (n * M_PI); + + gain /= fmax; // normalize + + for (int i = 0; i < ntaps; i++) + taps[i] *= gain; + + + return taps; +} + + vector gr_firdes::high_pass (double gain, - double sampling_freq, - double cutoff_freq, // Hz center of transition band - double transition_width, // Hz width of transition band - win_type window_type, - double beta) // used only with Kaiser + double sampling_freq, + double cutoff_freq, // Hz center of transition band + double transition_width, // Hz width of transition band + win_type window_type, + double beta) // used only with Kaiser { sanity_check_1f (sampling_freq, cutoff_freq, transition_width); int ntaps = compute_ntaps (sampling_freq, transition_width, - window_type, beta); + window_type, beta); // construct the truncated ideal impulse response times the window function @@ -129,15 +225,15 @@ gr_firdes::high_pass (double gain, taps[n + M] = -sin (n * fwT0) / (n * M_PI) * w[n + M]; } } - + // find the factor to normalize the gain, fmax. // For high-pass, gain @ fs/2 freq = 1.0 - + double fmax = taps[0 + M]; for (int n = 1; n <= M; n++) fmax += 2 * taps[n + M] * cos (n * M_PI); - gain /= fmax; // normalize + gain /= fmax; // normalize for (int i = 0; i < ntaps; i++) taps[i] *= gain; @@ -146,9 +242,57 @@ gr_firdes::high_pass (double gain, } // -// === Band Pass === +// === Band Pass === // +vector +gr_firdes::band_pass_2 (double gain, + double sampling_freq, + double low_cutoff_freq, // Hz center of transition band + double high_cutoff_freq, // Hz center of transition band + double transition_width, // Hz width of transition band + double attenuation_dB, // attenuation dB + win_type window_type, + double beta) // used only with Kaiser +{ + sanity_check_2f (sampling_freq, + low_cutoff_freq, + high_cutoff_freq, transition_width); + + int ntaps = compute_ntaps_windes (sampling_freq, transition_width, + attenuation_dB); + + vector taps(ntaps); + vector w = window (window_type, ntaps, beta); + + int M = (ntaps - 1) / 2; + double fwT0 = 2 * M_PI * low_cutoff_freq / sampling_freq; + double fwT1 = 2 * M_PI * high_cutoff_freq / sampling_freq; + + for (int n = -M; n <= M; n++){ + if (n == 0) + taps[n + M] = (fwT1 - fwT0) / M_PI * w[n + M]; + else { + taps[n + M] = (sin (n * fwT1) - sin (n * fwT0)) / (n * M_PI) * w[n + M]; + } + } + + // find the factor to normalize the gain, fmax. + // For band-pass, gain @ center freq = 1.0 + + double fmax = taps[0 + M]; + for (int n = 1; n <= M; n++) + fmax += 2 * taps[n + M] * cos (n * (fwT0 + fwT1) * 0.5); + + gain /= fmax; // normalize + + for (int i = 0; i < ntaps; i++) + taps[i] *= gain; + + return taps; +} + + vector gr_firdes::band_pass (double gain, double sampling_freq, @@ -201,6 +345,47 @@ gr_firdes::band_pass (double gain, // === Complex Band Pass === // +vector +gr_firdes::complex_band_pass_2 (double gain, + double sampling_freq, + double low_cutoff_freq, // Hz center of transition band + double high_cutoff_freq, // Hz center of transition band + double transition_width, // Hz width of transition band + double attenuation_dB, // attenuation dB + win_type window_type, + double beta) // used only with Kaiser +{ + sanity_check_2f_c (sampling_freq, + low_cutoff_freq, + high_cutoff_freq, transition_width); + + int ntaps = compute_ntaps_windes (sampling_freq, transition_width, + attenuation_dB); + + + + vector taps(ntaps); + vector lptaps(ntaps); + vector w = window (window_type, ntaps, beta); + + lptaps = low_pass_2(gain,sampling_freq,(high_cutoff_freq - low_cutoff_freq)/2,transition_width,attenuation_dB,window_type,beta); + + gr_complex *optr = &taps[0]; + float *iptr = &lptaps[0]; + float freq = M_PI * (high_cutoff_freq + low_cutoff_freq)/sampling_freq; + float phase=0; + if (lptaps.size() & 01) { + phase = - freq * ( lptaps.size() >> 1 ); + } else phase = - freq/2.0 * ((1 + 2*lptaps.size()) >> 1); + for(unsigned int i=0;i gr_firdes::complex_band_pass (double gain, double sampling_freq, @@ -240,11 +425,59 @@ gr_firdes::complex_band_pass (double gain, return taps; } - // // === Band Reject === // +vector +gr_firdes::band_reject_2 (double gain, + double sampling_freq, + double low_cutoff_freq, // Hz center of transition band + double high_cutoff_freq, // Hz center of transition band + double transition_width, // Hz width of transition band + double attenuation_dB, // attenuation dB + win_type window_type, + double beta) // used only with Kaiser +{ + sanity_check_2f (sampling_freq, + low_cutoff_freq, + high_cutoff_freq, transition_width); + + int ntaps = compute_ntaps_windes (sampling_freq, transition_width, + attenuation_dB); + + // construct the truncated ideal impulse response times the window function + + vector taps(ntaps); + vector w = window (window_type, ntaps, beta); + + int M = (ntaps - 1) / 2; + double fwT0 = 2 * M_PI * low_cutoff_freq / sampling_freq; + double fwT1 = 2 * M_PI * high_cutoff_freq / sampling_freq; + + for (int n = -M; n <= M; n++){ + if (n == 0) + taps[n + M] = 1.0 + ((fwT0 - fwT1) / M_PI * w[n + M]); + else { + taps[n + M] = (sin (n * fwT0) - sin (n * fwT1)) / (n * M_PI) * w[n + M]; + } + } + + // find the factor to normalize the gain, fmax. + // For band-reject, gain @ zero freq = 1.0 + + double fmax = taps[0 + M]; + for (int n = 1; n <= M; n++) + fmax += 2 * taps[n + M]; + + gain /= fmax; // normalize + + for (int i = 0; i < ntaps; i++) + taps[i] *= gain; + + return taps; +} + vector gr_firdes::band_reject (double gain, double sampling_freq, @@ -427,6 +660,19 @@ static const float width_factor[5] = { // indexed by win_type 10.0 // WIN_KAISER }; +int +gr_firdes::compute_ntaps_windes(double sampling_freq, + double transition_width, // this is frequency, not relative frequency + double attenuation_dB) +{ + // Based on formula from Multirate Signal Processing for + // Communications Systems, fredric j harris + int ntaps = (int)(attenuation_dB*sampling_freq/(22.0*transition_width)); + if ((ntaps & 1) == 0) // if even... + ntaps++; // ...make odd + return ntaps; +} + int gr_firdes::compute_ntaps (double sampling_freq, double transition_width, diff --git a/gnuradio-core/src/lib/general/gr_firdes.h b/gnuradio-core/src/lib/general/gr_firdes.h index f920398ad..bb62dad8f 100644 --- a/gnuradio-core/src/lib/general/gr_firdes.h +++ b/gnuradio-core/src/lib/general/gr_firdes.h @@ -68,6 +68,15 @@ class gr_firdes { win_type window = WIN_HAMMING, double beta = 6.76); // used only with Kaiser + static std::vector + low_pass_2 (double gain, + double sampling_freq, + double cutoff_freq, // Hz beginning transition band + double transition_width, // Hz width of transition band + double attenuation_dB, // attenuation dB + win_type window = WIN_HAMMING, + double beta = 6.76); // used only with Kaiser + /*! * \brief use "window method" to design a high-pass FIR filter * @@ -90,6 +99,15 @@ class gr_firdes { win_type window = WIN_HAMMING, double beta = 6.76); // used only with Kaiser + static std::vector + high_pass_2 (double gain, + double sampling_freq, + double cutoff_freq, // Hz center of transition band + double transition_width, // Hz width of transition band + double attenuation_dB, // attenuation dB + win_type window = WIN_HAMMING, + double beta = 6.76); // used only with Kaiser + /*! * \brief use "window method" to design a band-pass FIR filter * @@ -114,6 +132,15 @@ class gr_firdes { win_type window = WIN_HAMMING, double beta = 6.76); // used only with Kaiser + static std::vector + band_pass_2 (double gain, + double sampling_freq, + double low_cutoff_freq, // Hz beginning transition band + double high_cutoff_freq, // Hz beginning transition band + double transition_width, // Hz width of transition band + double attenuation_dB, // attenuation dB + win_type window = WIN_HAMMING, + double beta = 6.76); // used only with Kaiser /*! * \brief use "window method" to design a complex band-pass FIR filter @@ -140,6 +167,16 @@ class gr_firdes { win_type window = WIN_HAMMING, double beta = 6.76); // used only with Kaiser + static std::vector + complex_band_pass_2 (double gain, + double sampling_freq, + double low_cutoff_freq, // Hz beginning transition band + double high_cutoff_freq, // Hz beginning transition band + double transition_width, // Hz width of transition band + double attenuation_dB, // attenuation dB + win_type window = WIN_HAMMING, + double beta = 6.76); // used only with Kaiser + /*! * \brief use "window method" to design a band-reject FIR filter @@ -166,6 +203,16 @@ class gr_firdes { win_type window = WIN_HAMMING, double beta = 6.76); // used only with Kaiser + static std::vector + band_reject_2 (double gain, + double sampling_freq, + double low_cutoff_freq, // Hz beginning transition band + double high_cutoff_freq, // Hz beginning transition band + double transition_width, // Hz width of transition band + double attenuation_dB, // attenuation dB + win_type window = WIN_HAMMING, + double beta = 6.76); // used only with Kaiser + /*!\brief design a Hilbert Transform Filter * * \p ntaps: Number of taps, must be odd @@ -221,6 +268,10 @@ private: static int compute_ntaps (double sampling_freq, double transition_width, win_type window_type, double beta); + + static int compute_ntaps_windes (double sampling_freq, + double transition_width, + double attenuation_dB); }; #endif diff --git a/gnuradio-core/src/lib/general/qa_gr_firdes.cc b/gnuradio-core/src/lib/general/qa_gr_firdes.cc index 9184c8c92..5a41f118a 100644 --- a/gnuradio-core/src/lib/general/qa_gr_firdes.cc +++ b/gnuradio-core/src/lib/general/qa_gr_firdes.cc @@ -275,6 +275,211 @@ const static float t3_exp[107] = { }; +const static float t4_exp[] = { // low pass + 0.001059958362, +0.0002263929928, +-0.001277606934, +-0.0009675776237, + 0.001592264394, + 0.00243603508, +-0.001451682881, +-0.004769335967, +5.281541594e-18, + 0.007567512803, + 0.003658855334, +-0.009761494584, + -0.01011830103, + 0.009636915289, + 0.0193619132, +-0.004935568199, + -0.03060629964, +-0.007267376408, + 0.04236677289, + 0.03197422624, + -0.05274848267, + -0.0850463286, + 0.05989059806, + 0.31065014, + 0.4370569289, + 0.31065014, + 0.05989059806, + -0.0850463286, + -0.05274848267, + 0.03197422624, + 0.04236677289, +-0.007267376408, + -0.03060629964, +-0.004935568199, + 0.0193619132, + 0.009636915289, + -0.01011830103, +-0.009761494584, + 0.003658855334, + 0.007567512803, +5.281541594e-18, +-0.004769335967, +-0.001451682881, + 0.00243603508, + 0.001592264394, +-0.0009675776237, +-0.001277606934, +0.0002263929928, + 0.001059958362, +}; + + +const static float t5_exp[] = { //high pass +-0.001062123571, +-0.0002268554381, + 0.001280216733, + 0.000969554123, +-0.001595516922, +-0.002441011136, + 0.001454648213, + 0.004779078532, +-5.292330097e-18, +-0.007582970895, + -0.00366632943, + 0.009781434201, + 0.01013896987, +-0.009656600654, + -0.01940146461, + 0.004945650231, + 0.03066881932, + 0.00728222169, + -0.04245331511, + -0.03203954175, + 0.05285623297, + 0.08522006124, + -0.06001294032, + -0.3112847209, + 0.5630782247, + -0.3112847209, + -0.06001294032, + 0.08522006124, + 0.05285623297, + -0.03203954175, + -0.04245331511, + 0.00728222169, + 0.03066881932, + 0.004945650231, + -0.01940146461, +-0.009656600654, + 0.01013896987, + 0.009781434201, + -0.00366632943, +-0.007582970895, +-5.292330097e-18, + 0.004779078532, + 0.001454648213, +-0.002441011136, +-0.001595516922, + 0.000969554123, + 0.001280216733, +-0.0002268554381, +-0.001062123571, +}; + +const static float t6_exp[] = { // bandpass +0.0002809273137, +-0.001047327649, +7.936541806e-05, +-0.0004270860809, +0.0007595835486, +0.0008966081077, +-0.0004236323002, +0.0002423936094, +-0.002212299034, +0.0004807534278, +0.0002620361629, + 0.001443728455, + 0.002229931997, +-0.002720607212, +5.731141573e-05, +-0.004297634587, + 0.001878833398, + 0.003217151389, + 0.001357055153, + 0.003965090029, +-0.008576190099, +-0.0003257228818, +-0.004805727862, + 0.004721920472, + 0.01007549558, +-0.002688719891, + 0.004467967432, + -0.01837076992, +5.119658377e-17, + 0.001125075156, + 0.008071650751, + 0.02113764361, + -0.01602453552, + 0.001618095324, + -0.03004053794, + 0.003163811285, + 0.0219683405, + 0.007950295694, + 0.03682873398, + -0.05142467469, + -0.00794606097, + -0.03965795785, + 0.01544955093, + 0.09681399167, + -0.01610304788, + 0.08297294378, + -0.2811714709, + -0.1094062924, + 0.5275565982, + -0.1094062924, + -0.2811714709, + 0.08297294378, + -0.01610304788, + 0.09681399167, + 0.01544955093, + -0.03965795785, + -0.00794606097, + -0.05142467469, + 0.03682873398, + 0.007950295694, + 0.0219683405, + 0.003163811285, + -0.03004053794, + 0.001618095324, + -0.01602453552, + 0.02113764361, + 0.008071650751, + 0.001125075156, +5.119658377e-17, + -0.01837076992, + 0.004467967432, +-0.002688719891, + 0.01007549558, + 0.004721920472, +-0.004805727862, +-0.0003257228818, +-0.008576190099, + 0.003965090029, + 0.001357055153, + 0.003217151389, + 0.001878833398, +-0.004297634587, +5.731141573e-05, +-0.002720607212, + 0.002229931997, + 0.001443728455, +0.0002620361629, +0.0004807534278, +-0.002212299034, +0.0002423936094, +-0.0004236323002, +0.0008966081077, +0.0007595835486, +-0.0004270860809, +7.936541806e-05, +-0.001047327649, +0.0002809273137, +}; + void qa_gr_firdes::t1 () { @@ -340,5 +545,72 @@ qa_gr_firdes::t3 () void qa_gr_firdes::t4 () +{ + vector taps = + gr_firdes::low_pass_2 ( 1.0, + 8000, + 1750, + 500, + 66, + gr_firdes::WIN_HAMMING); + + // std::cout << "ntaps: " << taps.size () << std::endl; + // print_taps (std::cout, taps); + + CPPUNIT_ASSERT_EQUAL (NELEM (t4_exp), taps.size ()); + for (unsigned int i = 0; i < taps.size (); i++) + CPPUNIT_ASSERT_DOUBLES_EQUAL (t4_exp[i], taps[i], 1e-9); + + + check_symmetry (taps); +} + +void +qa_gr_firdes::t5 () +{ + vector taps = + gr_firdes::high_pass_2 ( 1.0, + 8000, + 1750, + 500, + 66, + gr_firdes::WIN_HAMMING); + + // std::cout << "ntaps: " << taps.size () << std::endl; + // print_taps (std::cout, taps); + + CPPUNIT_ASSERT_EQUAL (NELEM (t5_exp), taps.size ()); + + for (unsigned int i = 0; i < taps.size (); i++) + CPPUNIT_ASSERT_DOUBLES_EQUAL (t5_exp[i], taps[i], 1e-9); + + check_symmetry (taps); +} + +void +qa_gr_firdes::t6 () +{ + vector taps = + gr_firdes::band_pass_2 ( 1.0, + 20e6, + 5.75e6 - (5.28e6/2), + 5.75e6 + (5.28e6/2), + 0.62e6, + 66, + gr_firdes::WIN_HAMMING); + + // std::cout << "ntaps: " << taps.size () << std::endl; + // print_taps (std::cout, taps); + + CPPUNIT_ASSERT_EQUAL (NELEM (t6_exp), taps.size ()); + + for (unsigned int i = 0; i < taps.size (); i++) + CPPUNIT_ASSERT_DOUBLES_EQUAL (t6_exp[i], taps[i], 1e-7); + + check_symmetry (taps); +} + +void +qa_gr_firdes::t7 () { } diff --git a/gnuradio-core/src/lib/general/qa_gr_firdes.h b/gnuradio-core/src/lib/general/qa_gr_firdes.h index cb57fd740..d8566454a 100644 --- a/gnuradio-core/src/lib/general/qa_gr_firdes.h +++ b/gnuradio-core/src/lib/general/qa_gr_firdes.h @@ -32,6 +32,9 @@ class qa_gr_firdes : public CppUnit::TestCase { CPPUNIT_TEST (t2); CPPUNIT_TEST (t3); CPPUNIT_TEST (t4); + CPPUNIT_TEST (t5); + CPPUNIT_TEST (t6); + CPPUNIT_TEST (t7); CPPUNIT_TEST_SUITE_END (); private: @@ -39,6 +42,9 @@ class qa_gr_firdes : public CppUnit::TestCase { void t2 (); void t3 (); void t4 (); + void t5 (); + void t6 (); + void t7 (); }; -- cgit