<commit_before>/*
 * Copyright (c) 2007, 2008, 2009, 2010, 2011 <NAME>
 * Copyright (c) 2007, 2008, 2009, 2010, 2011 Virginia Polytechnic
 *                                      Institute & State University
 *
 * This file is part of liquid.
 *
 * liquid is free software: you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation, either version 3 of the License, or
 * (at your option) any later version.
 *
 * liquid is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU General Public License for more details.
 *
 * You should have received a copy of the GNU General Public License
 * along with liquid.  If not, see <http://www.gnu.org/licenses/>.
 */

//
// Finite impulse response filter design
//
// References:
//  [Herrmann:1973] <NAME>, <NAME>, and <NAME>,
//      "Practical design rules for optimum finite impulse response
//      lowpass digital filters," Bell Syst. Tech. Journal, vol. 52,
//      pp. 769--99, July-Aug. 1973
//  [Vaidyanathan:1993] <NAME>., "Multirate Systems and
//      Filter Banks," 1993, Prentice Hall, Section 3.2.1

#include <stdio.h>
#include <stdlib.h>
#include <math.h>

#include "liquid.internal.h"

// method to estimate required filter length
#define ESTIMATE_REQ_FILTER_LEN_METHOD_KAISER   (0)
#define ESTIMATE_REQ_FILTER_LEN_METHOD_HERRMANN (1)

#define ESTIMATE_REQ_FILTER_LEN_METHOD ESTIMATE_REQ_FILTER_LEN_METHOD_KAISER

// esimate required filter length given transition bandwidth and
// stop-band attenuation
//  _df     :   transition bandwidth (0 < _df < 0.5)
//  _As     :   stopband suppression level [dB] (_As > 0)
unsigned int estimate_req_filter_len(float _df,
                                     float _As)
{
    if (_df > 0.5f || _df <= 0.0f) {
        fprintf(stderr,"error: estimate_req_filter_len(), invalid bandwidth : %f\n", _df);
        exit(1);
    } else if (_As <= 0.0f) {
        fprintf(stderr,"error: estimate_req_filter_len(), invalid stopband level : %f\n", _As);
        exit(1);
    }

    // compute filter length estimate
#if ESTIMATE_REQ_FILTER_LEN_METHOD == ESTIMATE_REQ_FILTER_LEN_METHOD_KAISER
    // use Kaiser's estimate
    unsigned int h_len = (unsigned int) estimate_req_filter_len_Kaiser(_df,_As);
#elif ESTIMATE_REQ_FILTER_LEN_METHOD == ESTIMATE_REQ_FILTER_LEN_METHOD_HERRMANN
    // use Herrmann's estimate
    unsigned int h_len = (unsigned int) estimate_req_filter_len_Herrmann(_df,_As);
#else
#   error "invalid required filter length estimation method"
#endif
    
    return h_len;
}

// estimate filter stop-band attenuation given
//  _df     :   transition bandwidth (0 < _b < 0.5)
//  _N      :   filter length
float estimate_req_filter_As(float _df,
                             unsigned int _N)
{
    // run search for stop-band attenuation which gives these results
    float As0   = 0.01f;    // lower bound
    float As1   = 200.0f;   // upper bound

    float As_hat = 0.0f;    // stop-band attenuation estimate
    float N_hat = 0.0f;     // filter length estimate

    // perform simple bisection search
    unsigned int num_iterations = 20;
    unsigned int i;
    for (i=0; i<num_iterations; i++) {
        // bisect limits
        As_hat = 0.5f*(As1 + As0);
#if ESTIMATE_REQ_FILTER_LEN_METHOD == ESTIMATE_REQ_FILTER_LEN_METHOD_KAISER
        N_hat = estimate_req_filter_len_Kaiser(_df, As_hat);
#elif ESTIMATE_REQ_FILTER_LEN_METHOD == ESTIMATE_REQ_FILTER_LEN_METHOD_HERRMANN
        N_hat = estimate_req_filter_len_Herrmann(_df, As_hat);
#else
#       error "invalid required filter length estimation method"
#endif

        //printf("range[%8.2f, %8.2f] As-hat=%8.2fdB, N=%8.2f (target: %3u taps)\n",
        //        As0, As1, As_hat, N_hat, _N);

        // update limits
        if (N_hat < (float)_N) {
            As0 = As_hat;
        } else {
            As1 = As_hat;
        }
    }
    return As_hat;
}

// estimate filter transition bandwidth given
//  _As     :   stop-band attenuation [dB] (As > 0)
//  _N      :   filter length
float estimate_req_filter_df(float _As,
                             unsigned int _N)
{
    // run search for stop-band attenuation which gives these results
    float df0   = 1e-3f;    // lower bound
    float df1   = 0.499f;   // upper bound

    float df_hat = 0.0f;    // stop-band attenuation estimate
    float N_hat = 0.0f;     // filter length estimate

    // perform simple bisection search
    unsigned int num_iterations = 20;
    unsigned int i;
    for (i=0; i<num_iterations; i++) {
        // bisect limits
        df_hat = 0.5f*(df1 + df0);
#if ESTIMATE_REQ_FILTER_LEN_METHOD == ESTIMATE_REQ_FILTER_LEN_METHOD_KAISER
        N_hat = estimate_req_filter_len_Kaiser(df_hat, _As);
#elif ESTIMATE_REQ_FILTER_LEN_METHOD == ESTIMATE_REQ_FILTER_LEN_METHOD_HERRMANN
        N_hat = estimate_req_filter_len_Herrmann(df_hat, _As);
#else
#       error "invalid required filter length estimation method"
#endif

        printf("range[%8.5f, %8.5f] df-hat=%8.5fdB, N=%8.2f (target: %3u taps)\n",
                df0, df1, df_hat, N_hat, _N);

        // update limits
        if (N_hat < (float)_N) {
            df1 = df_hat;
        } else {
            df0 = df_hat;
        }
    }
    return df_hat;

}


// esimate required filter length given transition bandwidth and
// stop-band attenuation (algorithm from [Vaidyanathan:1993])
//  _df     :   transition bandwidth (0 < _df < 0.5)
//  _As     :   stop-band attenuation [dB] (As > 0)
float estimate_req_filter_len_Kaiser(float _df,
                                     float _As)
{
    if (_df > 0.5f || _df <= 0.0f) {
        fprintf(stderr,"error: estimate_req_filter_len_Kaiser(), invalid bandwidth : %f\n", _df);
        exit(1);
    } else if (_As <= 0.0f) {
        fprintf(stderr,"error: estimate_req_filter_len(), invalid stopband level : %f\n", _As);
        exit(1);
    }

    // compute filter length estimate
    return (_As - 7.95f)/(14.26f*_df);
}


// esimate required filter length given transition bandwidth and
// stop-band attenuation (algorithm from [Herrmann:1973])
//  _df     :   transition bandwidth (0 < _df < 0.5)
//  _As     :   stop-band attenuation [dB] (As > 0)
float estimate_req_filter_len_Herrmann(float _df,
                                       float _As)
{
    if (_df > 0.5f || _df <= 0.0f) {
        fprintf(stderr,"error: estimate_req_filter_len_Herrmann(), invalid bandwidth : %f\n", _df);
        exit(1);
    } else if (_As <= 0.0f) {
        fprintf(stderr,"error: estimate_req_filter_len(), invalid stopband level : %f\n", _As);
        exit(1);
    }

    // compute delta_1, delta_2
    float d1, d2;
    d1 = d2 = powf(10.0, -_As/10.0);

    // compute log of delta_1, delta_2
    float t1 = log10f(d1);
    float t2 = log10f(d2);

    // compute D_infinity(delta_1, delta_2)
    float Dinf = (0.005309f*t1*t1 + 0.07114f*t1 - 0.4761f)*t2 -
                 (0.002660f*t1*t1 + 0.59410f*t1 + 0.4278f);

    // compute f(delta_1, delta_2)
    float f = 11.012f + 0.51244f*(t1-t2);

    // compute filter length estimate
    float h_len = (Dinf - f*_df*_df) / _df + 1.0f;
    
    return h_len;
}<commit_msg>firdes|kaiser window: fixing bug in Herrmann's filter length computation<commit_after>/*
 * Copyright (c) 2007, 2008, 2009, 2010, 2011 <NAME>
 * Copyright (c) 2007, 2008, 2009, 2010, 2011 Virginia Polytechnic
 *                                      Institute & State University
 *
 * This file is part of liquid.
 *
 * liquid is free software: you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation, either version 3 of the License, or
 * (at your option) any later version.
 *
 * liquid is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU General Public License for more details.
 *
 * You should have received a copy of the GNU General Public License
 * along with liquid.  If not, see <http://www.gnu.org/licenses/>.
 */

//
// Finite impulse response filter design
//
// References:
//  [Herrmann:1973] <NAME>, <NAME>, and <NAME>,
//      "Practical design rules for optimum finite impulse response
//      lowpass digital filters," Bell Syst. Tech. Journal, vol. 52,
//      pp. 769--99, July-Aug. 1973
//  [Vaidyanathan:1993] <NAME>., "Multirate Systems and
//      Filter Banks," 1993, Prentice Hall, Section 3.2.1

#include <stdio.h>
#include <stdlib.h>
#include <math.h>

#include "liquid.internal.h"

// method to estimate required filter length
#define ESTIMATE_REQ_FILTER_LEN_METHOD_KAISER   (0)
#define ESTIMATE_REQ_FILTER_LEN_METHOD_HERRMANN (1)

// select filter estimate method
#define ESTIMATE_REQ_FILTER_LEN_METHOD          (0)

// esimate required filter length given transition bandwidth and
// stop-band attenuation
//  _df     :   transition bandwidth (0 < _df < 0.5)
//  _As     :   stopband suppression level [dB] (_As > 0)
unsigned int estimate_req_filter_len(float _df,
                                     float _As)
{
    if (_df > 0.5f || _df <= 0.0f) {
        fprintf(stderr,"error: estimate_req_filter_len(), invalid bandwidth : %f\n", _df);
        exit(1);
    } else if (_As <= 0.0f) {
        fprintf(stderr,"error: estimate_req_filter_len(), invalid stopband level : %f\n", _As);
        exit(1);
    }

    // compute filter length estimate
#if ESTIMATE_REQ_FILTER_LEN_METHOD == ESTIMATE_REQ_FILTER_LEN_METHOD_KAISER
    // use Kaiser's estimate
    unsigned int h_len = (unsigned int) estimate_req_filter_len_Kaiser(_df,_As);
#elif ESTIMATE_REQ_FILTER_LEN_METHOD == ESTIMATE_REQ_FILTER_LEN_METHOD_HERRMANN
    // use Herrmann's estimate
    unsigned int h_len = (unsigned int) estimate_req_filter_len_Herrmann(_df,_As);
#else
#   error "invalid required filter length estimation method"
#endif
    
    return h_len;
}

// estimate filter stop-band attenuation given
//  _df     :   transition bandwidth (0 < _b < 0.5)
//  _N      :   filter length
float estimate_req_filter_As(float _df,
                             unsigned int _N)
{
    // run search for stop-band attenuation which gives these results
    float As0   = 0.01f;    // lower bound
    float As1   = 200.0f;   // upper bound

    float As_hat = 0.0f;    // stop-band attenuation estimate
    float N_hat = 0.0f;     // filter length estimate

    // perform simple bisection search
    unsigned int num_iterations = 20;
    unsigned int i;
    for (i=0; i<num_iterations; i++) {
        // bisect limits
        As_hat = 0.5f*(As1 + As0);
#if ESTIMATE_REQ_FILTER_LEN_METHOD == ESTIMATE_REQ_FILTER_LEN_METHOD_KAISER
        N_hat = estimate_req_filter_len_Kaiser(_df, As_hat);
#elif ESTIMATE_REQ_FILTER_LEN_METHOD == ESTIMATE_REQ_FILTER_LEN_METHOD_HERRMANN
        N_hat = estimate_req_filter_len_Herrmann(_df, As_hat);
#else
#       error "invalid required filter length estimation method"
#endif

        //printf("range[%8.2f, %8.2f] As-hat=%8.2fdB, N=%8.2f (target: %3u taps)\n",
        //        As0, As1, As_hat, N_hat, _N);

        // update limits
        if (N_hat < (float)_N) {
            As0 = As_hat;
        } else {
            As1 = As_hat;
        }
    }
    return As_hat;
}

// estimate filter transition bandwidth given
//  _As     :   stop-band attenuation [dB] (As > 0)
//  _N      :   filter length
float estimate_req_filter_df(float _As,
                             unsigned int _N)
{
    // run search for stop-band attenuation which gives these results
    float df0   = 1e-3f;    // lower bound
    float df1   = 0.499f;   // upper bound

    float df_hat = 0.0f;    // stop-band attenuation estimate
    float N_hat = 0.0f;     // filter length estimate

    // perform simple bisection search
    unsigned int num_iterations = 20;
    unsigned int i;
    for (i=0; i<num_iterations; i++) {
        // bisect limits
        df_hat = 0.5f*(df1 + df0);
#if ESTIMATE_REQ_FILTER_LEN_METHOD == ESTIMATE_REQ_FILTER_LEN_METHOD_KAISER
        N_hat = estimate_req_filter_len_Kaiser(df_hat, _As);
#elif ESTIMATE_REQ_FILTER_LEN_METHOD == ESTIMATE_REQ_FILTER_LEN_METHOD_HERRMANN
        N_hat = estimate_req_filter_len_Herrmann(df_hat, _As);
#else
#       error "invalid required filter length estimation method"
#endif

        //printf("range[%8.5f, %8.5f] df-hat=%8.5fdB, N=%8.2f (target: %3u taps)\n",
        //        df0, df1, df_hat, N_hat, _N);

        // update limits
        if (N_hat < (float)_N) {
            df1 = df_hat;
        } else {
            df0 = df_hat;
        }
    }
    return df_hat;

}


// esimate required filter length given transition bandwidth and
// stop-band attenuation (algorithm from [Vaidyanathan:1993])
//  _df     :   transition bandwidth (0 < _df < 0.5)
//  _As     :   stop-band attenuation [dB] (As > 0)
float estimate_req_filter_len_Kaiser(float _df,
                                     float _As)
{
    if (_df > 0.5f || _df <= 0.0f) {
        fprintf(stderr,"error: estimate_req_filter_len_Kaiser(), invalid bandwidth : %f\n", _df);
        exit(1);
    } else if (_As <= 0.0f) {
        fprintf(stderr,"error: estimate_req_filter_len(), invalid stopband level : %f\n", _As);
        exit(1);
    }

    // compute filter length estimate
    return (_As - 7.95f)/(14.26f*_df);
}


// esimate required filter length given transition bandwidth and
// stop-band attenuation (algorithm from [Herrmann:1973])
//  _df     :   transition bandwidth (0 < _df < 0.5)
//  _As     :   stop-band attenuation [dB] (As > 0)
float estimate_req_filter_len_Herrmann(float _df,
                                       float _As)
{
    if (_df > 0.5f || _df <= 0.0f) {
        fprintf(stderr,"error: estimate_req_filter_len_Herrmann(), invalid bandwidth : %f\n", _df);
        exit(1);
    } else if (_As <= 0.0f) {
        fprintf(stderr,"error: estimate_req_filter_len(), invalid stopband level : %f\n", _As);
        exit(1);
    }

    // compute delta_1, delta_2
    float d1, d2;
    d1 = d2 = powf(10.0, -_As/20.0);

    // compute log of delta_1, delta_2
    float t1 = log10f(d1);
    float t2 = log10f(d2);

    // compute D_infinity(delta_1, delta_2)
    float Dinf = (0.005309f*t1*t1 + 0.07114f*t1 - 0.4761f)*t2 -
                 (0.002660f*t1*t1 + 0.59410f*t1 + 0.4278f);

    // compute f(delta_1, delta_2)
    float f = 11.012f + 0.51244f*(t1-t2);

    // compute filter length estimate
    float h_len = (Dinf - f*_df*_df) / _df + 1.0f;
    
    return h_len;
}
{'fraction_non_alphanumeric': 0.10951893551688843, 'fraction_numerical': 0.0373251450017059, 'mean_word_length': 2.549527730685396, 'pattern_counts': {'":': 0, '<': 57, '<?xml version=': 0, '>': 35, 'https://': 0, 'lorem ipsum': 0, 'www.': 2, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 8, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '3708004', 'n_tokens_mistral': 5329, 'n_tokens_neox': 4863, 'n_words': 1372}