Фильтрация с помощью БПФ

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

При увеличении порядка фильтра (увеличении длины линии задержки) число умножений растёт.

Когда порядок фильтра превышает 50, становится выгоднее использовать расчёт отклика фильтра не прямым способом, а через умножение спектров в частотной области с последующим обратным БПФ.

Известны 2 метода секционированной свёртки: метод перекрытия с накоплением и метод перекрытия с суммированием.

Выходной буфер имеет размер (длина БПФ - (длина импульсной характеристики фильтра -1)).

С помощью данных методов можно обрабатывать сразу 2 канала. Для этого потребуется использовать комплексное БПФ.

Для этого буфер для импульсной характеристики длиной L должен иметь в 2 раза больший размер. Отбрасываться будут (L - 1) отсчётов для каждого канала, то есть всего (L - 1) * 2 отсчётов.

Буфер для импульсной характеристики перед комплексным БПФ

i[0], 0, i[1], 0, i[2], 0, ..., i[imp.length - 1], 0, 0, 0 ...

Буфер, содержащий входные данные, готовится в виде:

l[0], r[0], l[1], r[1], ..., l[n - 1], r[n - 1]

Выходной буфер комплексного БПФ будет иметь размер (1 << power) - 2 * (L - 1).

Выходные данные в нём расположены так же, как и во входной буфере:

l[0], r[0], l[1], r[1], ..., l[(1 << power) - 2 * (L - 1) - 1], r[(1 << power) - N * (L - 1) - 1]

Расчёт импульсной характеристики фильтра можно использовать из раздела КИХ фильтры.

Естественно, можно использовать и АЧХ фильтра.

Реализация метода перекрытия с накоплением

/**

* фильтрация с помощью БПФ

* алгоритм взят из

* Романюк Ю.А. Основы цифровой обработки сигналов.

* 5.4 Секционированная свёртка, Метод перекрытия с накоплением.

*/

public class FFTfilter1 {

 /**

  * импульсная характеристика после преобразования в частотную область

  */

 private double[] m_response;

 /**

  * буфер для накопления нужного числа входных данных

  */

 private short[] m_buffer;

 /**

  * указатель для сохранения входных данных в буфер

  */

 private int m_s;

 /**

  * размер буфера, степень числа 2

  */

 private byte m_power;

 /**

  * длина импульсной характеристики фильтра - 1

  */

 private int m_imp_length1;

 //

 public FFTfilter1() {}

/**

* инициализация параметров

* @param imp - импульсная характеристика фильтра

*/

 public void init(double[] imp) {

         // преобразуем импульсную характеристику в частотную область

         // подбираем размер буфера для свёртки.

         // (L - 1) отсчётов будут отбрасываться,

         // длину буфера выбираем минимум 32768.

         int i;

         m_power = 0;

         for( i = 1; i < imp.length; m_power++ ) i = i << 1;

         // если степень меньше 15, возрастают искажения

         if( ++m_power < 15 ) m_power = 15;

         // создаём буфер для преобразования и дополняем его нулями

         m_response = new double[1 << m_power];// уже обнулён

         for( i = 0; i < imp.length; i++ ) m_response[i] = imp[i];

         //for( ; i < (1 << m_power); i++ ) m_response[i] = 0.0;

         // переводим в частотную область

         FFT.realFastFourierTransform( m_response, m_power, false );

         // вспомогательный буфер для хранения входных данных

         m_buffer = new short[1 << m_power];

         m_imp_length1 = imp.length - 1;

         m_s = m_imp_length1;// компенсируем начальное смещение влево

         // обнуляем кусок буфера, в java он уже обнулён

         //for( i = 0; i < m_s; i++ ) m_buffer[i] = 0;

 }

/**

* фильтрация

* выходной сигнал имеет задержку на L/2 сэмплов

* @param sample - отсчёт звука

* @return результат фильтрации, если не null

*/

 public short[] proc(short sample) {

         m_buffer[m_s++] = sample;

         if( m_s < m_buffer.length )

                 return null;

         // переносим данные во временный буфер

         double[] buf = new double[m_buffer.length];

         for( int i = 0; i < m_buffer.length; i++ )

                 buf[i] = (double)m_buffer[i];

         // сохраняем последние (L - 1) значений для следующего преобразования

         for( int i = m_buffer.length - m_imp_length1, m = 0; i < m_buffer.length; i++ )

                 m_buffer[m++] = m_buffer[i];

         m_s = m_imp_length1;

         // выполняем БПФ

         FFT.realFastFourierTransform( buf, m_power, false );

         // перемножаем частотные характеристики

         /* постоянная составляющая

         buf[0] = buf[0] * m_response[0];

         buf[1] = buf[1] * m_response[1];

         for( int m = 2; m < buf.length; m += 2 )

         */

         for( int m = 0; m < buf.length; m += 2 ) {

                 // от способа перемножения зависит,

                 // какая часть данных будет отбрасываться

                 double a = buf[m];

                 double b = buf[m + 1];

                 // вариант для отбрасывания последних (L-1) отсчётов

                 buf[m]     = a * m_response[m] + b * m_response[m + 1];

                 buf[m + 1] = b * m_response[m] - a * m_response[m + 1];

                 /* вариант для отбрасывания первых (L-1) отсчётов

                 buf[m]     = a * m_response[m] - b * m_response[m + 1];

                 buf[m + 1] = b * m_response[m] + a * m_response[m + 1];

                 */

         }

         // преобразуем обратно во временную область

         FFT.realFastFourierTransform( buf, m_power, true );

         // переписываем часть данных в выходной буфер

         short[] out = new short[m_buffer.length - m_imp_length1];

         // отбрасываем последние (L-1) отсчётов

         for( int i = 0; i < out.length; i++ ) {

                 if( buf[i] > 32767 ) out[i] = 32767;

                 else if( buf[i] < -32768 ) out[i] = -32768;

                 else out[i] = (short)buf[i];

         }

         /* или отбрасываем первые (L-1) отсчётов

         for( int i = 0; i < out.length; i++ ) {

                 if( buf[i + m_imp_length1] > 32767 ) out[i] = 32767;

                 else if( buf[i + m_imp_length1] < -32768 ) out[i] = -32768;

                 else out[i] = (short)buf[i + m_imp_length1];

         }

         */

         return out;

 }

}

Реализация метода перекрытия с суммированием

/**

* фильтрация с помощью БПФ

* алгоритм взят из

* Романюк Ю.А. Основы цифровой обработки сигналов.

* 5.4 Секционированная свёртка, Метод перекрытия с суммированием.

*/

public class FFTfilter2 {

 /**

  * импульсная характеристика после преобразования в частотную область

  */

 private double[] m_response;

 /**

  * буфер для накопления нужного числа входных данных

  */

 private short[] m_buffer;

 /**

  * указатель для сохранения входных данных в буфер

  */

 private int m_s;

 /**

  * размер буфера, степень числа 2

  */

 private byte m_power;

 /**

  * длина импульсной характеристики фильтра - 1

  */

 private int m_imp_length1;

 /**

  * буфер для хранения выходных данных для последующего суммирования

  */

 private double[] m_out;

 //

 public FFTfilter2() {}

/**

* инициализация параметров

* @param imp - импульсная характеристика фильтра

*/

 public void init(double[] imp) {

         // преобразуем импульсную характеристику в частотную область

         // подбираем размер буфера для свёртки.

         // (L - 1) отсчётов будут использоваться для суммирования,

         // длину буфера выбираем минимум 32768.

         int i;

         m_power = 0;

         for( i = 1; i < imp.length; m_power++ ) i = i << 1;

         // если степень меньше 15, возрастают искажения

         if( ++m_power < 15 ) m_power = 15;

         // создаём буфер для преобразования и дополняем его нулями

         m_response = new double[1 << m_power];// уже обнулён

         for( i = 0; i < imp.length; i++ ) m_response[i] = imp[i];

         //for( ; i < (1 << m_power); i++ ) m_response[i] = 0.0;

         // переводим в частотную область

         FFT.realFastFourierTransform( m_response, m_power, false );

         // вспомогательный буфер для хранения входных данных

         m_imp_length1 = imp.length - 1;

         m_buffer = new short[(1 << m_power) - m_imp_length1];

         m_s = 0;

         m_out = new double[m_imp_length1];// заполнен нулями

 }

/**

* фильтрация

* выходной сигнал имеет задержку,

* зависящую от длины импульсной характеристики фильтра

* @param sample - отсчёт звука

* @return результат фильтрации, если не null

*/

 public short[] proc(short sample) {

         m_buffer[m_s++] = sample;

         if( m_s < m_buffer.length )

                 return null;

         m_s = 0;

         // переносим данные во временный буфер

         double[] buf = new double[1 << m_power];

         int i = 0;

         do{ buf[i] = (double)m_buffer[i]; } while( ++i < m_buffer.length );

         // обнуляем последние (L-1) значений

         do{ buf[i] = 0; } while( ++i < buf.length );// в java это не требуется

         // выполняем БПФ

         FFT.realFastFourierTransform( buf, m_power, false );

         // перемножаем частотные характеристики

         /* постоянная составляющая

         buf[0] = buf[0] * m_response[0];

         buf[1] = buf[1] * m_response[1];

         for( i = 2; i < buf.length; i += 2 )

         */

         for( i = 0; i < buf.length; i += 2 ) {

                 double a = buf[i];

                 double b = buf[i + 1];

                 buf[i]     = a * m_response[i] - b * m_response[i + 1];

                 buf[i + 1] = a * m_response[i + 1] + b * m_response[i];

         }

         // преобразуем обратно во временную область

         FFT.realFastFourierTransform( buf, m_power, true );

         // переписываем часть данных в выходной буфер

         // суммируем последние (L-1) отсчётов старого буфера

         // и первые (L-1) отсчётов нового буфера

         short[] out = new short[buf.length - m_imp_length1];

         for( i = 0; i < m_imp_length1; i++ ) {

                 double v = m_out[i] + buf[i];

                 if( v > 32767 ) out[i] = 32767;

                 else if( v < -32768 ) out[i] = -32768;

                 else out[i] = (short)v;

         }

         for( ; i < out.length; i++ ) {

                 if( buf[i] > 32767 ) out[i] = 32767;

                 else if( buf[i] < -32768 ) out[i] = -32768;

                 else out[i] = (short)buf[i];

         }

         // сохраняем последние (L-1) значений для следующего суммирования

         for( ; i < buf.length; i++ ) { m_out[i - out.length] = buf[i]; }

         return out;

 }

}