Изменение частоты дискретизации

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

Ошибка складывается из двух составляющих:

- ошибка позиционирования на временной оси

- ошибка аппроксимации

Первая составляющая легко уменьшается путем увеличения точности расчета. Для того, чтобы добиться уровня влияния ниже -90 дБ необходимо иметь точность не ниже 16-ти двоичных разрядов. Это условие обычно легко выполнимо.

Вторая составляющая зависит от выбранного способа аппроксимации. Самый простой способ - линейная аппроксимация. Он даёт искажения примерно на уровне -60 дБ для частоты 1кГц при передискретизации с 44100 Гц на 48000 Гц. Интерполяция по Лагранжу даёт искажения на уровне -90дБ.

Analog Devices в своих преобразователях частоты дискретизации использует следующий способ:

- поиск положения на временной оси осуществляется с точностью 20 двоичных разрядов

- аппроксимация (восстановление огибающей) осуществляется с помощью фильтрации 64-х точечным КИХ-фильтром.

Эти операции также соответствуют следующему порядку действий: сначала в 1000000 раз увеличили частоту дискретизации и восстановили огибающую, затем из полученных отсчётов выбрали ближайший к интересующей точке времени. Искажения при таких параметрах находятся на уровне -120 дБ. Прямая реализация такого алгоритма невозможна, так как частота дискретизации получается несколько десятков гигагерц и требуется производить много вычислений, результат которых будет просто выбрасываться, поэтому используются различные ухищрения.

Кроме прямого назначения, этот же алгоритм может применяться для изменения времени звучания звукового фрагмента, аналогично притормаживанию-ускорению пластинки или плёнки при проигрывании.

Алгоритм линейной аппроксимации для изменения частоты дискретизации применён для получения значений сэмплированного голосового источника для любой частоты дискретизации. См. Голосовой сэмплированный источник.

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

Спектр 1кГц после линейной интерполяции

Спектр 1кГц после целочисленной интерполяции по Лагранжу с коэффициентом 3

Спектр 7кГц после линейной интерполяции

Спектр 7кГц после целочисленной интерполяции по Лагранжу с коэффициентом 3

Спектр 7кГц после интерполяции КИХ фильтром с 10 умножениями на отсчёт и коэффициентом передискретизации 1024

Спектр 7кГц после интерполяции КИХ фильтром с 20 умножениями на отсчёт, коэффициентом передискретизации 32 и интерполяцией коэффициентов фильтра

Спектр 18кГц после интерполяции КИХ фильтром с 10 умножениями на отсчёт и коэффициентом передискретизации 16384

Спектр 18кГц после интерполяции КИХ фильтром с 40 умножениями на отсчёт и коэффициентом передискретизации 16384

Спектр 18кГц после интерполяции КИХ фильтром с 80 умножениями на отсчёт, коэффициентом передискретизации 32 и интерполяцией значений фильтра

Спектр 18кГц после интерполяции КИХ фильтром со 100 умножениями на отсчёт и коэффициентом передискретизации 16384

Смотри также

Передискретизация

Secret Rabbit Code

http://www.mega-nerd.com/SRC/

Реализация для линейной интерполяции

/**

* Преобразователь частоты дискретизации

* с линейной интерполяцией.

* Выходной сигнал имеет задержку на 1 сэмпл

*

*/

public class SRCLinear {

 private byte m_accuracy;

 private int m_mask;

 private int m_dTf;

 private int m_dT;

 private int m_Tf;

 private int m_T;

 private int m_inT;

 private short m_previous;

 //

 public SRCLinear() {

         init( 1, 1, (byte)16 );

 }

 public boolean init(int inSampleRate, int outSampleRate, byte accuracy) {

         m_mask = (1 << accuracy) - 1;

         long quotient = ((long)inSampleRate << accuracy) / (long)outSampleRate;

         m_dTf = (int)(quotient & m_mask);

         m_dT = (int)(quotient >> accuracy);

         m_accuracy = accuracy;

         m_Tf = 0;

         m_T = 0;

         m_inT = 0;

         m_previous = 0;

         return (quotient > 0);

 }

 /**

  * @param in: входной буфер

  * @param out: выходной буфер, должен иметь достаточный размер

  * @return количество сэмплов в выходном буфере

  */

 public int process(short[] in, short[] out) {

         int outCount = 0;

         int s = 0;

         do {

                 while( m_inT <= m_T ) {

                         m_previous = in[s++];

                         m_inT++;

                         if( s >= in.length ) {

                                 if( m_inT < m_T ) {

                                         m_T = m_T - m_inT;

                                         m_inT = 0;

                                 } else {

                                         m_inT = m_inT - m_T;

                                         m_T = 0;

                                 }

                                 return outCount;

                         }

                 }

                 out[outCount++] = (short)((int)m_previous + (int)(((long)m_Tf * (long)((int)in[s] - (int)m_previous)) >> m_accuracy));

                 long t = (long)m_Tf + (long)m_dTf;

                 m_T = m_T + m_dT + (int)(t >> m_accuracy);

                 m_Tf = (int)(t & m_mask);

         }

         while( outCount < out.length );

         if( m_inT < m_T ) {

                 m_T = m_T - m_inT;

                 m_inT = 0;

         } else {

                 m_inT = m_inT - m_T;

                 m_T = 0;

         }

         return outCount;

 }

}

Реализация для интерполяции по Лагранжу

/**

* Преобразователь частоты дискретизации

* с интерполяцией по Лагранжу.

* Выходной сигнал имеет задержку на (N+1)/2 сэмпл.

*

*/

public class SRCLagrange {

 private double m_dT;

 private double m_T;

 private int m_inT;

 private short[] m_pF;

 private byte m_wrPos;

 private double[] m_pLI;

 private byte m_interpolation;

 private byte m_Hinterpolation;

 //

 public SRCLagrange() {

         init( 1, 1, (byte)1 );

 }

 public boolean init(int inSampleRate, int outSampleRate, byte interpolation)

 {

         if( (interpolation & 1) == 0 ) return false;

         m_interpolation = interpolation;

         m_Hinterpolation = (byte)((interpolation - 1) >> 1);

         m_pF = new short[m_interpolation];

         m_pLI = new double[m_interpolation];

         m_wrPos = m_interpolation;

         m_dT = (double)inSampleRate / (double)outSampleRate;

         m_T = 0.0;

         m_inT = 0;

         return (m_dT > 0.0);

 }

 /**

  * @param in: входной буфер

  * @param out: выходной буфер, должен иметь достаточный размер

  * @return количество сэмплов в выходном буфере

  */

 public int process(short[] in, short[] out) {

         int outCount = 0;

         int s = 0;

         do {

                 while( m_inT <= (int)m_T ) {

                         m_inT++;

                         m_pF[m_wrPos++] = in[s++];

                         if( m_wrPos > m_interpolation ) m_wrPos = 0;

                         if( s >= in.length ) {

                                 if( m_inT < m_T ) {

                                         m_T = m_T - (double)m_inT;

                                         m_inT = 0;

                                 } else {

                                         m_inT = m_inT - (int)m_T;

                                         m_T = m_T - (double)(int)m_T;

                                 }

                                 return outCount;

                         }

                 }

                 double D = (double)m_Hinterpolation + m_T - (double)(int)m_T;

                 for( byte n = 0; n <= m_interpolation; n++ ) {

                         m_pLI[n] = 1.0;

                         for( byte k = 0; k <= m_interpolation; k++ )

                                 if( n != k )

                                         m_pLI[n] = m_pLI[n] * (D - (double)k) / (double)(n - k);

                 }

                 double fout = 0.0;

                 for( byte k = 0, index = m_wrPos; ; ) {

                         fout = fout + m_pLI[k++] * m_pF[index];

                         if( k > m_interpolation ) break;

                         if( ++index > m_interpolation ) index = 0;

                 }

                 // ограничение амплитуды

                 int iout = (int)fout;

                 if( iout > 32767 ) iout = 32767;

                 else if( iout < -32768 ) iout = -32768;

                 out[outCount++] = (short)iout;

                 //

                 m_T = m_T + m_dT;

         }

         while( outCount < out.length );

         if( m_inT < m_T ) {

                 m_T = m_T - (double)m_inT;

                 m_inT = 0;

         } else {

                 m_inT = m_inT - (int)m_T;

                 m_T = m_T - (double)(int)m_T;

         }

         return outCount;

 }

}

Реализация для интерполяции c помощью КИХ-фильтра

/**

* Преобразователь частоты дискретизации

* с помощью КИХ-фильтра.

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

* длины фильтра.

*

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

* oversampling можно добавить вычисление методом линейной интерполяции

* отсутствующих значений КИХ фильтра.

*

*/

public class SRCFIR {

 private byte m_mul_count;

 private int m_oversampling;

 private int m_fir_length1;

 private double[] m_fir;

 private short[] m_buff;

 private double m_dT;

 private double m_T;

 private int m_inT;

 /**

  * описание параметров смотрите в {@link #init(int, int, int, byte)}

  */

 public SRCFIR(int inSampleRate, int outSampleRate, int oversampling, byte mul)

 {

         init( inSampleRate, outSampleRate, oversampling, mul );

 }

 /**

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

  * @param inSampleRate - частота дискретизации входного сигнала

  * @param outSampleRate - частота дискретизации выходного сигнала

  * @param oversampling - частота дискретизации фильтра,

  *     определяет точность по оси времени.

  *     если не используется интерполяция коэффициентов фильтра,

  *     стоит использовать значение не менее 1024.

  *     в аппаратных ПЧД используется 1048576.

  * таблица коэффициентов фильтра имеет размер oversampling * mul.

  * произведение oversampling * mul должно быть чётным для

  * расчёта фильтра чётного порядка, что позволяет хранить

  * только половину коэффициентов и съэкономить память.

  * @param mul - желаемое число умножений при фильтрации.

  *     определяет загрузку процессора.

  * @return true, если успешна, иначе false

  */

 public boolean init(int inSampleRate, int outSampleRate, int oversampling, byte mul)

 {

         if( (mul & oversampling & 1) != 0 )

                 return false;

         m_oversampling = oversampling;

         m_mul_count = mul;

         m_fir_length1 = m_oversampling * m_mul_count;

         m_fir = new double[m_fir_length1 >> 1];

         // расчитываем фильтр для половины частоты дискретизации

         // так как рабочая частота повышена в oversampling раз,

         // частота фильтра должна быть в oversampling раз меньше.

         // f = 1, fd = 2

         // w = 2 * pi * f / fd / oversampling

         double w = Math.PI / m_oversampling;

         double sum = 0.0;

         // для симметричной характеристики

         double c = (double)m_fir.length - 0.5;

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

                 double d = (double)i - c;

                 m_fir[i] = Math.sin(w * d) / d * FIR.getWindow( i, m_fir_length1, FIR.BLACKMAN );

                 sum += m_fir[i] + m_fir[i];

         }

         m_fir_length1--;

         // нормализация и масштабирование

         sum /= oversampling;

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

                 m_fir[i] /= sum;

         m_buff = new short[m_mul_count];

         for( int i = m_mul_count; i-- > 0; ) m_buff[i] = 0;

         m_dT = (double)inSampleRate / (double)outSampleRate;

         m_T = 0.0;

         m_inT = 0;

         return true;

 }

 /**

  * @param in: входной буфер

  * @param out: выходной буфер, должен иметь достаточный размер

  * @return количество сэмплов в выходном буфере

  */

 public int process(short[] in, short[] out) {

         int outCount = 0;

         int s = 0;

         do {

                 while( m_inT <= m_T ) {

                         m_inT++;

                         // буфер для фильтрации, линия задержки

                         for( int i = m_mul_count; --i > 0; )

                                 m_buff[i] = m_buff[i - 1];

                         m_buff[0] = in[s++];

                         if( s >= in.length ) {

                                 if( m_inT < m_T ) {

                                         m_T = m_T - (double)m_inT;

                                         m_inT = 0;

                                 } else {

                                         m_inT = m_inT - (int)m_T;

                                         m_T = m_T - (double)(int)m_T;

                                 }

                                 return outCount;

                         }

                 }

                 // ближайшее начальное смещение по временной оси

                 int shift = (int)(0.5 + m_oversampling * (m_T - (double)(int)m_T));

                 // для увеличения точности при небольшом значении oversampling

                 // можно дополнительно расчитать методом интерполяции

                 // значение коэффицента фильтра в требуемой точке.

                 // фильтрация

                 double fout = 0.0;

                 int k, i;

                 // первая половина КИХ-фильтра

                 for( k = shift, i = 0; k < m_fir.length; k += m_oversampling )

                         fout += m_fir[k] * m_buff[i++];

                 // вторая половина КИХ-фильтра

                 for( k = m_fir_length1 - k; k >= 0; k -= m_oversampling )

                         fout += m_fir[k] * m_buff[i++];

                 // ограничение амплитуды

                 if( fout >= 0 ) {

                         fout += 0.5;

                         if( fout > 32767.0 ) fout = 32767.0;

                 } else {

                         fout -= 0.5;

                         if( fout < -32768.0 ) fout = -32768.0;

                 }

                 out[outCount++] = (short)fout;

                 //

                 m_T = m_T + m_dT;

         }// на всякий случай проверяем не переполнился ли выходной буфер

         while( outCount < out.length );

         if( m_inT < m_T ) {

                 m_T = m_T - (double)m_inT;

                 m_inT = 0;

         } else {

                 m_inT = m_inT - (int)m_T;

                 m_T = m_T - (double)(int)m_T;

         }

         return outCount;

 }

}