Различные выходы при применении фильтра к аудиосигналу

Question

Goodevening, я пишу простое программное обеспечение на C, которое читает аудиофайл wav с использованием библиотеки libsndfile (http://www.mega-nerd.com/libsndfile/api.html), чем образцы переходят к функции процесса, где я применяю фильтр к сигнал (фильтр низких частот Баттерворта второго порядка, применяется с транспонированной прямой формой II). После этого я пишу результат в новый wav-файл. Если вместо фильтра я применяю простые операции (например, для muliply samples для константы), он отлично работает, но когда я применяю фильтр, он создает много шума. Я попытался напечатать значения образцов после фильтра и до того, как они будут записаны в новом файле, и я получил те же значения, что и у Matlab (где результат, который я получаю, является идеальным), но они отличаются от значений Я получил, если я прочитал вывод, написанный библиотекой.

static void processAudio (double *buffer, int length) 
{ 
    //arrays a and b are the coefficients 
    double a[] = { 1, 
        -1.799096409484668, 
        0.817512403384758}; 


    double b[] = { 0.004603998475022, 
        0.009207996950045, 
        0.004603998475022}; 

    //arrays a_ and b_ are the coefficients normalized 
    double a_[] = {1,a[1]/b[0],a[2]/b[0]}; 
    double b_[] = {1,b[1]/b[0],b[2]/b[0]}; 
    double gain = b[0]/a[0]; 

    double reg[] ={0,0}; //memory registers 

    for(int i = 0; i<length; i++) 
    { 
     if(i%2==0) //just left channel is changed 
     { 
      //TRANSPOSED DIRECT FORM II 
      double input = *(buffer+i); 
      double output =(input + reg[0])*gain; 

      reg[0] = reg[1]+b_[1]*input-a_[1]*output; 
      reg[1] = b_[2]*input - a_[2]*output; 

      *(buffer+i) = output; 

     } 
    } 
} 



//that function is called in the main method inside this cycle 
int main(void) 
{ 
    ... 

    while ((framesRead = sf_read_double(inputFile, buffer, BUFFER_LENGTH))) 
    { 

      processAudio(buffer, framesRead); 

      sf_write_double(outputFile, buffer, framesRead); 
    } 

    ... 
} 

Если напечатать результаты только после фильтра я получаю:

0.00000000000000 
0.00000000000000 
0.00000000000000 
0.00000000000000 
0.00000000000000 
0.00000000000000 
0.00000000000000 
0.00000000000000 
0.00000000000000 
0.00000000000000 
0.00000000000000 
0.00000000000000 
0.00000000000000 
0.00000000000000 
0.00000000000000 
-0.00000014050288 
-0.00000025277823 
0.00000050310774 
0.00000209530885 
0.00000420138311 
0.00000725078112 
0.00001115570320 
0.00001554761090 
0.00002053775981 
0.00002550357112 
0.00002796948679 
0.00000727086200 
-0.00006401853354 
... 

Если я прочитаю выходной файл результаты

0 
0 
0 
0 
0 
0 
0 
0 
0 
0 
0 
0 
0 
0 
0 
0 
0 
0 
0 
0 
0 
0 
3.05175781250000e-05 
3.05175781250000e-05 
3.05175781250000e-05 
3.05175781250000e-05 
0 
-6.10351562500000e-05 
-0.000366210937500000 
-0.00183105468750000 
-0.00601196289062500 
-0.0140075683593750 
-0.0261840820312500 
-0.0421142578125000 
-0.0610656738281250 
-0.0820617675781250 
-0.104156494140625 
-0.126342773437500 
-0.147766113281250 
-0.167663574218750 
-0.185363769531250 
-0.200378417968750 
-0.212280273437500 
... 

Они очень разные. Я не знаю, что происходит. Если у вас есть идея, дайте мне знать. Заранее спасибо!

Answer 1

Вы получаете разные выходы, потому что не сохраняете состояния фильтра. На каждой итерации цикла кадров вы возвращаете состояние фильтра в начальное состояние {0.0, 0.0}. Кроме того, не нормализуйте коэффициенты, так как MATLAB делает это для вас. Просто возьмите коэффициенты из MATLAB и примените его к фильтру.

Я сформирует-структуру, чтобы сохранить состояние фильтра, например:

struct biquad { 
    double b0, b1, b2, a1, a2; // Note: MATLAB will always normalize a0 to 1.0, so no need to process that. 
    double r0, r1; 
}; 

Затем я хотел бы добавить функцию как biquad_process, что-то вроде этого:

void process_biquad(struct biquad *self, double *buffer, int length) 
{ 
    int i; 

    for (i = 0; i < length; i++) { 
     double x = buffer[i]; 
     double y = (x * self->b0) + self->r0; 

     self->r0 = self->r1 + (x * self->b1) - (self->a1 * y); 
     self->r1 = (self->b2 * x) - (self->a2 * y); 

     buffer[i] = y; 
    } 
} 

Тогда, в вашей функции main() вы можете загрузить аудио в буфере и обработать каждый фрагмент следующим образом:

int main(void) 
{ 
    struct biquad *bq = calloc(1, sizeof(struct biquad)); 

    // Take coefficients from MATLAB and put them here 
    bq->b0 = ... 
    bq->b1 = ... 

    // Set initial state to 0.0 
    bq->r0 = bq->r1 = 0.0; 
    ... 

    while ((framesRead = sf_read_double(inputFile, buffer, BUFFER_LENGTH))) 
    { 

      process_biquad(bq, buffer, framesRead); 

      sf_write_double(outputFile, buffer, framesRead); 
    } 

    ... 
} 

Надеюсь, это поможет.