EarFilter.cpp

Go to the documentation of this file.

#include "EarFilter.h"
 
//=====================================================================================
krotos::EarFilter::EarFilter() {}
 
//=====================================================================================
void krotos::EarFilter::configure(float sampleRate)
{
    m_sampleRate = sampleRate;
 
    // extend to sampleRate / 2
    if (m_sampleRate / 2.0f > m_maxFreq)
    {
        m_freqs.push_back(m_sampleRate / 2.0f); // extend to sr / 2
        m_gains.push_back(m_gains.back());
    }
 
    // second order low pass coefficients
    m_a = getLowPassFilterCoefficients(m_fc / m_sampleRate, m_q);
 
    // convert to highpass
    m_b = convertToHighPass(m_a);
 
    // create the filter object and filter the signal
    // create the filter object and filter the signal
    // the only way it does not crash
    ReferenceCountedObjectPtr<CoefficientsIIR> filterCoeff =
            new CoefficientsIIR(m_b[0], m_b[1], m_b[2], m_a[0], m_a[1], m_a[2]);
    m_filter.coefficients = filterCoeff;
 
    // transfer function of filter applied so far
    calculateTransferFunction(m_b, m_a, 512);
 
    // square every element of mag response
    std::transform(m_magResponse.begin(), m_magResponse.end(), m_magResponse.begin(),
                   [](float value) { return value * value; });
    // transfer function that remains to apply
    std::vector<float> gainsInterp = LinearInterpolator::interpolate(m_freqs, m_gains, m_normalisedFreqs);
    // divide by transfer function of first filter
    m_gain = elementWiseDivision(gainsInterp, m_magResponse);
 
    // apply gain condition for elements below m_freqs(2)
    float conditionValue = m_freqs[1];
 
    for (size_t i = 0; i < m_gain.size(); ++i)
    {
        if (m_normalisedFreqs[i] < conditionValue)
        {
            m_gain[i] = 1.0f;
        }
    }
 
    m_freqResponse = linspace(0.0f, 1.0f, static_cast<int>(m_gain.size()));
 
    // numerator
    m_firCoeffb = fir2(N, m_freqResponse, m_gain);
}
 
//=====================================================================================
std::vector<float> krotos::EarFilter::processSignal(std::vector<float>& inputSignal)
{
    jassert(m_sampleRate != -1.0f);
 
    const int numSamples = static_cast<int>(inputSignal.size());
 
    // filter the signal twice in cascade
    std::vector<float> filteredSignal1;
    std::vector<float> filteredSignal2;
    applyIIRFilter(inputSignal, m_filter, filteredSignal1);
    applyIIRFilter(filteredSignal1, m_filter, filteredSignal2);
 
    // Apply to signal
    std::vector<float> filteredSignal3;
    applyFIRFilter(filteredSignal2, m_firCoeffb, filteredSignal3);
 
    return filteredSignal3;
}
 
std::vector<float> krotos::EarFilter::getLowPassFilterCoefficients(float normalizedFrequency, float qualityFactor)
{
    jassert(qualityFactor != 0.0f);
 
    float rho = std::exp(-float(M_PI) * normalizedFrequency / qualityFactor);
    float theta =
            2.0f * float(M_PI) * normalizedFrequency * std::sqrt(1.0f - 1.0f / (4.0f * std::pow(qualityFactor, 2.0f)));
 
    std::vector<float> coeff(static_cast<int>(rho), 1.0f);
 
    float a1 = -2.0f * rho * cos(theta);
    float a2 = std::pow(rho, 2.0f);
 
    coeff.push_back(1.0f);
    coeff.push_back(a1);
    coeff.push_back(a2);
 
    return coeff;
}
 
std::vector<float> krotos::EarFilter::convertToHighPass(const std::vector<float>& lowPassCoefficients)
{
    float sum = 0.0f;
    for (int i = 0; i < lowPassCoefficients.size(); i++)
    {
        sum += lowPassCoefficients[i];
    }
 
    float b0 = sum - 1.0f;
    float b1 = -lowPassCoefficients[1];
    float b2 = -lowPassCoefficients[2];
 
    std::vector<float> coeff;
    coeff.push_back(b0);
    coeff.push_back(b1);
    coeff.push_back(b2);
 
    return coeff;
}
 
void krotos::EarFilter::applyIIRFilter(const std::vector<float>& inputSignal, dsp::IIR::Filter<float>& filter,
                                       std::vector<float>& filteredSignal)
{
    for (const auto& sample : inputSignal)
    {
        float filteredSample = filter.processSample(sample);
        filteredSignal.push_back(filteredSample);
    }
}
 
void krotos::EarFilter::applyFIRFilter(const std::vector<float>& inputSignal, const std::vector<float>& coefficients,
                                       std::vector<float>& filteredSignal)
{
    int numSamples = static_cast<int>(inputSignal.size());
    int numCoefficients = static_cast<int>(coefficients.size());
 
    for (int sample = 0; sample < numSamples; ++sample)
    {
        float filteredSample = 0.0f;
 
        for (int i = 0; i < numCoefficients; ++i)
        {
            int index = sample - i;
            if (index >= 0)
                filteredSample += coefficients[i] * inputSignal[index];
        }
 
        filteredSignal.push_back(filteredSample);
    }
}
 
void krotos::EarFilter::calculateTransferFunction(const std::vector<float> b, const std::vector<float> a, int numPoints)
{
    CoefficientsIIR coeff(b[0], b[1], b[2], a[0], a[1], a[2]);
 
    // evaluate the transfer function
    m_magResponse.resize(numPoints);
    m_normalisedFreqs.resize(numPoints);
 
    for (int i = 0; i < numPoints; ++i)
    {
        float normalizedFrequency = static_cast<float>(i) / numPoints;
        float magnitude = static_cast<float>(coeff.getMagnitudeForFrequency(
                static_cast<double>(normalizedFrequency * 0.5f * m_sampleRate), static_cast<double>(m_sampleRate)));
 
        // square magnitude values
        m_magResponse[i] = magnitude;
        // 0...pi to 0...Fs/2
        m_normalisedFreqs[i] = normalizedFrequency * m_sampleRate / 2.0f;
    }
}
 
std::vector<float> krotos::EarFilter::elementWiseDivision(const std::vector<float>& vector1,
                                                          const std::vector<float>& vector2)
{
    // check if input vectors are of equal size
    jassert(vector1.size() == vector2.size());
 
    std::vector<float> result(vector1.size()); // reserve memory for the result
 
    for (int i = 0; i < vector1.size(); ++i)
    {
        // perform element-wise division
        result[i] = (vector2[i] != 0.0f) ? (vector1[i] / vector2[i] + std::numeric_limits<float>::epsilon()) : 0.0f;
    }
 
    return result;
}
 
std::vector<float> krotos::EarFilter::fir2(int number, std::vector<float>& freqResponse, const std::vector<float>& gain)
{
    // check freq and gain should have same size
    jassert(freqResponse.size() == gain.size());
    // work with filter length
    number += 1;
    float npt = 512;
 
    if (number < 1024)
        npt = 512;
    else
        npt = std::pow(2.0f, std::ceilf(std::log2(static_cast<float>(number))));
 
    std::vector<float> wind(number);
    // create Window
    for (int n = 0; n < number; ++n)
    {
        wind[n] = 0.54f - 0.46f * std::cos((2.0f * float(M_PI) * n) / float(number - 1.0f));
    }
 
    // Ensure the window is symmetric
    for (int i = 0; i < number / 2; ++i)
    {
        wind[number - i - 1] = wind[i];
    }
 
    int lap = static_cast<int>(std::floor(npt / 25.0f));
 
    int nbrk = static_cast<int>(gain.size());
 
    const float eps = std::numeric_limits<float>::epsilon();
 
    jassert(std::abs(freqResponse[0]) <= eps || std::abs(freqResponse[nbrk - 1] - 1) <= eps);
 
    freqResponse[0] = 0.0f;
    freqResponse[nbrk - 1] = 1.0f;
 
    // interpolate breakpoints onto large grid
 
    int nint = nbrk - 1;
    std::vector<double> df(freqResponse.size());
    adjacent_difference(freqResponse.begin(), freqResponse.end(), df.begin());
    df.erase(df.begin());
 
    // jassert(std::any_of(df.begin(), df.end(), [](double val) { return val < 0;
    // }));
 
    if (number > 2 * static_cast<int>(npt))
    {
        jassertfalse;
    }
 
    // length of[dc 1 2 ... nyquist] frequencies.
    npt += 1.0f;
    std::vector<float> H;
    int nb = 1;
    int ne;
 
    for (int i = 0; i < nint; ++i)
    {
        if (df[i] == 0)
        {
            nb = static_cast<int>(std::ceil(nb - lap / 2.0f));
            ne = nb + lap;
        }
        else
        {
            ne = static_cast<int>(std::floor(freqResponse[i + 1] * npt));
        }
 
        if (nb < 0 || ne > npt)
        {
            jassertfalse;
        }
 
        for (int j = nb; j <= ne; ++j)
        {
            float inc = (j - nb) / static_cast<float>(ne - nb);
 
            if (nb == ne)
                inc = 0;
 
            H.push_back(inc * gain[i + 1] + (1.0f - inc) * gain[i]);
        }
 
        nb = ne + 1;
    }
 
    // Fourier time shift
    float dt = 0.5f * (number - 1);
 
    std::vector<std::complex<float>> rad(static_cast<int>(npt));
 
    for (int i = 0; i < npt; ++i)
    {
        float realPart = 0.0f;
        float imagPart = -dt * float(M_PI) * static_cast<float>(i) / static_cast<float>(npt - 1);
        ; // Since the imaginary part is 0 in this case
        rad[i] = std::complex<float>(realPart, imagPart);
    }
 
    std::vector<std::complex<float>> Hexp(H.size());
 
    for (int i = 0; i < H.size(); ++i)
    {
        std::complex<float> expValue = std::exp(rad[i]);
        Hexp[i] = H[i] * expValue;
    }
 
    // force last element to real
    Hexp.back().imag(0.0f);
 
    // Fourier transform of real series
    // reverse H excluding the last element
    std::vector<std::complex<float>> H_reversed(Hexp.rbegin() + 1, Hexp.rend());
    std::vector<std::complex<float>> H_conj;
    for (const auto& h : H_reversed)
    {
        // calculate the conjugate of each element
        H_conj.push_back(std::conj(h));
    }
 
    // concatenate H with its reversed conjugate
    Hexp.insert(Hexp.end(), H_conj.begin(), H_conj.end() - 1);
 
    // initialise FFT object
    juce::dsp::FFT fft(static_cast<int>(std::log2(static_cast<float>(Hexp.size()))));
 
    // perform the inverse FFT
    std::vector<std::complex<float>> Hexp1(Hexp.size());
    fft.perform(Hexp.data(), Hexp1.data(), true);
 
    // extract the real part of the result
    std::vector<float> ht(Hexp1.size());
 
    for (int i = 0; i < Hexp1.size(); ++i)
    {
        ht[i] = Hexp1[i].real();
    }
 
    // raw numerator
    std::vector<float> b(ht.begin(), ht.begin() + number);
 
    // apply Window
    for (int i = 0; i < number; i++)
    {
        b[i] *= wind[i];
    }
 
    return b;
}