Главная » Программирование звука » Как измерить одну частоту

0

Если у вас есть мензурка, полная морской воды, вы могли бы задаться вопросом: сколько хлорида калия в этой мензурке? Точно так же, если у нас есть оцифрованный  звук,  который  состоит  из  множества  различных  синусоид,  можно  было бы  поинтересоваться,  насколько  сильна  синусоида  частотой  1000  Гц  в  этом  звуке?  Дискретное  преобразование  Фурье  (ДПФ)  дает  возможность  ответить  на  этот вопрос.  Если  вам  интересна  только  одна  частота,  используйте  ДПФ  для  измерения  того,  насколько  сильна  эта  частота  как  часть  общей  звуковой  картины.  Повторяя  это  измерение  для  разных  частот,  вы  можете  получить  достаточно  полное представление о спектре звукового сигнала.

Хотя  ДПФ  опирается  на  достаточно  сложную  математику,  основная  идея  проста: если вы умножите две разных синусоиды, то получите разные результаты, которые зависят от отношения между частотами.

Например,  предположим,  что  вы  хотите  выяснить,  насколько  сильно  представлена  волна  частотой  1000  Гц  в  определенном  сигнале.  Сначала  умножаем  заданный сигнал на синусоиду с частотой 1000 Гц, а затем складываем вместе полученные отсчеты. Ha рис. 24.4 показан этот процесс для случая, когда исследуемый сигнал есть синусоида с частотой 1000 Гц с той же фазой. (Очевидно, что картина будет  меняться  в  зависимости  от  фазы   мы  подробно  рассмотрим  это  позднее.) Может это и неочевидно, но в нашем случае конечная сумма пропорциональна амплитуде исследуемого сигнала.

B противоположность этому, на рис. 24.5 показано, что произойдет, если иссле-

дуемый сигнал является синусоидой с отличной частотой. Здесь в произведении

умноженная  на  f/N.  Подобным  образом,  t   это  целочисленный  номер  отсчета. Кроме  того,  суммирование  дает  не  непосредственное  значение  амплитуды,  а  всего  лишь  число,  пропорциональное  амплитуде.  Несмотря  на  эти  оговорки,  результат все же заслуживает внимания.

Если  вы  повторите  эти  вычисления  для  различных  значений  f,  то  сможете  измерить  амплитуду  всех  частот  в вашем  сигнале.  Для  любого  целого  f меньшего  N вы   легко   определите   значение   Af,   представляющее   амплитуду   соответствующей частоты как долю от общего сигнала. Эти значения могут быть вычислены по той же формуле:

?1

Af   = ? st cos(2?tf  / N )

=0

Если   мы   знаем   значения   Аf,   мы   можем   восстановить   исходные   отсчеты. Вспомним,  что  Af     амплитуда  синусоиды  с  частотой  fS/N.  Для  восстановления сигнала  необходимо  сложить  все  значения  для  разных  частот.  Аналитически  это записывается так:

?1

st  = ? Af  cos(2?tf  / )

=0

Обратное преобразование  аналогично  прямому, о  котором мы  говорили выше. Если  вы  постоянно  работаете  с  преобразованиями  Фурье,  это  подобие  становится важным.

Учет фазы

Чтобы  осуществлять  точное  обратное  преобразование  Фурье,  помимо  амплитуды и частоты необходимо измерять фазу каждой частоты. Мы уже знаем, как измерять амплитуду каждой частоты. A фазу?

Ha  помощь  приходят  комплексные  числа.  Можно  изменить  описанный  ранее метод  вычислений  так,  что  он  будет  давать  двумерный  результат.  Простое  комплексное  число   это,  по  существу,  двумерное  значение,  поэтому  оно  одновременно представляет и амплитуду, и фазу.

При   таком   подходе   фазовая   часть   вычисляется   неявно.   Вместо   амплитуды и  фазы  вы  измеряете  две  амплитуды,  соответствующие  разным  фазам.  Одна  из этих фаз представляется косинусом (cos()), другая синусом (sin()).

Используя  комплексные  числа,  вы  можете  проводить  измерения  одновремен-

но,  умножая  синусную  часть  на  -i.  (Для  представления торые предпочитают букву j)

N ?1

? 1   я  использую  i,  неко-

Af   = ? (st cos(2?tf  / N ) ? ist sin(2?tf  / N ))

=0

Каждое   значение   Af    теперь   представляется   комплексным   числом;   действи-

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

фазами.  Вспомним,  что  любое  комплексное  число  a+ib  имеет  модуль,  равный

a 2  + b 2  ,  и  аргумент,  равный  tan-1(a/b).  Модуль  Af    это  амплитуда  одной  синусо-

идальной волны, фаза которой равна аргументу Af.

Для   компактности   обычно   заменяют   выражение   cos(2?tf/N)     isin(2?tf/N) на  эквивалентное  e-2?tf/N,  что  дает  вариант  записи  ДПФ,  обычно  используемый в книгах:

N ?1

Af   = ?

t =0

st e

? 2?itf  / N

He удивительно, что полное обратное ДПФ выглядит абсолютно так же. Теперь вы должны понимать, что при такой формулировке для каждой частоты происходит комбинирование  синусоиды  и  косинусоиды  для  получения  одной  волны  с  истинной амплитудой, частотой и фазой. Сложение этих волн дает вам исходный сигнал:

?1

s() = ?

=0

Af e

? 2?itf / N

Реализация ДПФ

Для  ДПФ  легко  разработать  программный  код.  Используя  библиотеку  complex из стандартной библиотеки C++, мы можем реализовать прямое ДПФ способом, приведенным в листинге 24.1, воспользовавшись тем, что polar(1.0, x) то же самое,

что и eix.

Листинг 24.1. Дискретное преобразование Фурье

void ForwardDft(complex <double> *samples,

int length, complex<double>*result) {

static const double

twoPi = 2*3.1415926535897932384626; for(int f = 0; f<length; f++) { result[f] = 0.0;

for(int t = 0; t < length; t++)

result[f] += samples[t] * polar(1.0,-twoPi*f*t/length);

}

}

Эта  функция  заманчиво  проста,  и  вам  следует  всегда  держать  ее  под  рукой, если  вы  часто  работаете  с  преобразованиями  Фурье.  Можно  не  сомневаться,  что благодаря  этой  простоте  она  работает  правильно.  Сравнивая  другие  программы преобразования Фурье с этой, вы можете проверять верность их работы.

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

Основная  проблема  заключается  в  том,  что  эта  формула  вычисляет  каждую точку  спектра  по  отдельности.  Ha  практике  вам  почти  всегда  требуется  обсчитать   сразу   весь   спектр,   что   позволяет   совместить   вычисления   для   различных частот.

Масштабирование

Если  вы  возьмете  набор  отсчетов,  примените  к  нему  прямое  ДПФ,  описанное ранее,  для  получения  спектра,  а  затем  используете  обратное  ДПФ  для  получения набора  отсчетов,  то  отсчеты  будут  иными  по  сравнению  с  исходными.  Результат в действительности будет в N раз больше. Чтобы скомпенсировать это, нужно гденибудь ввести поправочный множитель  Существует, по крайней мере, три способа сделать это.

Один  из  способов  получения  корректно  промасштабированного  результата  разделить  результаты  прямого  ДПФ  на  N.  Прямое  и  обратное  ДПФ  тогда  будут выглядеть так:

1 N ?1

A f   =

t =0

st e

? 2?itf  / N

N ?1

s() = ?

f =0

Af e

? 2?itf / N

Теперь,  если  вы  примените  прямое  ДПФ  к  набору  отсчетов,  а  затем  обратное

ДПФ к полученному спектру, результат, как уже упоминалось, будет нулевым.

Некоторые   предпочитают   делить   результаты   прямого   и   обратного   ДПФ

на    , что дает исключительно симметричные формулы:

A f   =

N ?1

s e ? 2?itf  / N

N  t =0

N ?1

s(t ) =

1   ? A e ?2?itf / N

N  f =0

Однако  обычно  программистов  заботит  время,  необходимое  для  выполнения прямого   преобразования,   поэтому   наиболее   распространенный   метод   корректировки заключается в делении результата обратного ДПФ на N:

N ?1

A   = ? s e ? 2?itf  / N

t =0

1 N ?1

f

N f =0

Хотя  первый  метод  и  претендует  на  некоторую  техническую  корректность, основной   практический   интерес   представляет   относительная   амплитуда   выходного  сигнала.  Любой  из  этих  трех  методов  масштабирования  дает  одни  и  те  же относительные результаты.

Источник: Кинтцель Т.  Руководство программиста по работе со звуком = A Programmer’s Guide to Sound: Пер. с англ. М.: ДМК Пресс, 2000. 432 с, ил. (Серия «Для программистов»).

По теме:

  • Комментарии