Реализация алгоритма Гёртцеля в C

Question

Я внедряю систему BFSK со скачкообразной перестройкой частоты на процессоре DSP. Некоторые участники форума предложили использовать алгоритм Гёртцеля для демодуляции скачкообразной перестройки частоты на определенных частотах. Я попытался реализовать алгоритм Goertzel в C. код выглядит следующим образом:

float goertzel(int numSamples,int TARGET_FREQUENCY,int SAMPLING_RATE, float* data)
{
    int     k,i;
    float   floatnumSamples;
    float   omega,sine,cosine,coeff,q0,q1,q2,result,real,imag;

    floatnumSamples = (float) numSamples;
    k = (int) (0.5 + ((floatnumSamples * TARGET_FREQUENCY) / SAMPLING_RATE));
    omega = (2.0 * M_PI * k) / floatnumSamples;
    sine = sin(omega);
    cosine = cos(omega);
    coeff = 2.0 * cosine;
    q0=0;
    q1=0;
    q2=0;

    for(i=0; i<numSamples; i++)
    {
        q0 = coeff * q1 - q2 + data[i];
        q2 = q1;
        q1 = q0;
    }
    real = (q1 - q2 * cosine);
    imag = (q2 * sine);
    result = sqrtf(real*real + imag*imag);
    return result;
}

Когда я использую функцию для вычисления результата на определенных частотах для данного набора данных, я не получаю правильные результаты. Однако, если я использую тот же набор данных и вычисляю результат goertzel с помощью функции MATLAB goertzel(), то я получаю результаты идеально. Я реализовал алгоритм, используя C, с помощью некоторых онлайн-учебников, которые я нашел через Интернет. Я просто хочу узнать ваше мнение, ребята, правильно ли функция реализует алгоритм Гёртцеля.

Answer 1

Если вы говорите, что реализация Matlab хороша, потому что ее результаты совпадают с результатами для этой частоты DFT или FFT ваших данных, то это, вероятно, потому, что реализация Matlab нормализует результаты с помощью коэффициента масштабирования, как это делается с FFT.

Измените код, чтобы учесть это и посмотреть, улучшит ли он ваши результаты. Обратите внимание, что я также изменил имена функций и результатов, чтобы отразить, что ваш goertzel вычисляет величину, а не полный комплексный результат, для ясности:

float goertzel_mag(int numSamples,int TARGET_FREQUENCY,int SAMPLING_RATE, float* data)
{
    int     k,i;
    float   floatnumSamples;
    float   omega,sine,cosine,coeff,q0,q1,q2,magnitude,real,imag;

    float   scalingFactor = numSamples / 2.0;

    floatnumSamples = (float) numSamples;
    k = (int) (0.5 + ((floatnumSamples * TARGET_FREQUENCY) / SAMPLING_RATE));
    omega = (2.0 * M_PI * k) / floatnumSamples;
    sine = sin(omega);
    cosine = cos(omega);
    coeff = 2.0 * cosine;
    q0=0;
    q1=0;
    q2=0;

    for(i=0; i<numSamples; i++)
    {
        q0 = coeff * q1 - q2 + data[i];
        q2 = q1;
        q1 = q0;
    }

    // calculate the real and imaginary results
    // scaling appropriately
    real = (q1 - q2 * cosine) / scalingFactor;
    imag = (q2 * sine) / scalingFactor;

    magnitude = sqrtf(real*real + imag*imag);
    return magnitude;
}

Answer 2

Часто вы можете просто использовать квадрат величины в своих вычислениях, например, для обнаружения тона.

Некоторые превосходные примеры Goertzels находятся в коде DSP Asterisk PBX, коде DSP Asterisk (dsp.c) и в библиотеке spandsp. Библиотека SPANDSP DSP.

Answer 3

Рассмотрим две входные выборочные формы волны:

1) синусоида с амплитудой А и частотой W

2) косинусная волна с одинаковой амплитудой и частотой A и W

Алгоритм Гёртцеля должен давать одинаковые результаты для двух упомянутых входных волновых форм, но предоставленный код приводит к различным возвращаемым значениям. Я думаю, что код должен быть пересмотрен следующим образом:

float goertzel_mag(int numSamples,int TARGET_FREQUENCY,int SAMPLING_RATE, float* data)
{
    int     k,i;
    float   floatnumSamples;
    float   omega,sine,cosine,coeff,q0,q1,q2,magnitude,real,imag;

    float   scalingFactor = numSamples / 2.0;

    floatnumSamples = (float) numSamples;
    k = (int) (0.5 + ((floatnumSamples * TARGET_FREQUENCY) / SAMPLING_RATE));
    omega = (2.0 * M_PI * k) / floatnumSamples;
    sine = sin(omega);
    cosine = cos(omega);
    coeff = 2.0 * cosine;
    q0=0;
    q1=0;
    q2=0;

    for(i=0; i<numSamples; i++)
    {
        q2 = q1;
        q1 = q0;
        q0 = coeff * q1 - q2 + data[i];
    }

    // calculate the real and imaginary results
    // scaling appropriately
    real = (q0 - q1 * cosine) / scalingFactor;
    imag = (-q1 * sine) / scalingFactor;

    magnitude = sqrtf(real*real + imag*imag);
    return magnitude;
}