Цель этой статьи не в том, чтобы дать строгое математическое определение преобразованию Фурье — это бесчисленное количество раз уже сделано другими авторами, а на примерах показать его «механический» смысл и пояснить почему оно работает.

Читателю понадобятся знания математики и физики на уровне 9 класса средней школы. Я предполагаю, что вы знаете теорему Пифагора, чему равен синус и косинус угла в прямоугольном треугольнике и что такое синусоида.

В конце статьи мы применим полученные знания для решения стандартной задачи декодирования телефонных номеров, сохранённых в аудиофайл в виде DTMF сигналов.

Безответственность

Категорически не рекомендую данную статью математикам и людям, предпочитающим формальную строгость и точность доказательств. Здесь нет никаких доказательств и никакой строгости. Вместо них — анимации, картинки, образные сравнения и наглядные примеры, понятные людям с неполным средним образованием. Сам автор едва владеет арифметикой и из всех цифр помнит только первые две.

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

Виртуальные приборы в нашей лаборатории

Для наших исследований мы будем использовать пару виртуальных самописцев. Их задача — принимать на вход некий сигнал, например, напряжение электрического тока в вольтах, и перемещать пишущий механизм от начального «нулевого» положения на расстояние, соответствующее значению измеряемого сигнала.

На рис. 1 показано устройство трехканального самописца, регистрирующего одновременно три разных величины.

Самописец состоит из пишущего узла, включающего в себя перо (О₁, О₂, О₃) с пропитанным краской пористым материалом внутри и электромеханического привода пера (П₁, П₂, П₃), перемещающего перо в зависимости от уровня подаваемого в привод электрического тока. Прямо под пишущим узлом обычно находится бумажная лента (М), которая с постоянной скоростью протягивается при помощи лентопротяжного механизма (Пр).

Во время работы самописца мы получаем на ленте развёртку входного сигнала во времени. Зная скорость перемещения ленты, можно посчитать, в какой момент времени подаваемый в самописец сигнал имел то или иное значение.

Другая разновидность самописцев работает не с лентой, а с бумажным кругом (рис. 2). Пишущий механизм у них аналогичен таковому в ленточных приборах, но в качестве носителя информации используется вращающийся бумажный круг. Это устройство мы тоже используем в наших исследованиях.

Опыт первый: Level Zero

Запустим наши виртуальные приборы и подадим на них нулевой сигнал. Пишущие перья установились в среднее положение, соответствующее нулю на входе. Скорость движения ленты — одно деление в секунду. Скорость вращения диска — один сектор в секунду, а всего секторов там 12.

Перья обоих приборов работают абсолютно синхронно, одинаково реагируя на входной сигнал. В данном случае они неподвижны. Из под пера ленточного самописца выходит по шесть точек за секунду: именно с такой скоростью наши приборы способны измерять и записывать уровень входного сигнала. Я специально добавил такую особенность к виртуальным самописцам чтобы сделать их поведение более «цифровым». Далее нам это потребуется.

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

Опыт второй: непрямая прямая

В прямой нет ничего особенного, вы знаете уравнение:

Попробуем изобразить на самописце прямую Я написал в формуле , потому что по оси абсцисс (горизонтальной) на ленточном самописце у нас, по сути, откладывается время по одному делению за одну секунду.

Таак... ну, левый самописец, можно сказать, рисует прямую... Но правый явно делает что-то не то! Не прямую, а какую-то...

гипножабу

Разница в получении этих двух изображений линейной функции заключается только в способе движения бумажки под перьями самописцев. Сами перья, очевидно, совершают абсолютно одинаковые перемещения.

Опыт третий: круглая синусоида

Подадим на самописцы синусоидальный сигнал:

Глядя на работу ленточного прибора, можно заметить, что полный период, равный радиан (360 градусов), синусоида проходит за 12 делений ленты или за те же самые 12 реальных секунд. Нетрудно вычислить частоту и угловую частоту:

Формула нашего сигнала:

В этом примере:

• амплитуда (измеряем в вертикальных делениях на ленте),

• период с,

• частота , Гц,

• угловая частота: ; радиан/с,

• начальная фаза равна нулю.

Итоговая формула

означает, что за время секунд аргумент функции изменяется на радиан (или 360 градусов), после чего, в силу периодичности синуса, значения начнут повторяться сначала.

С тем же самым периодом секунд вращается диск кругового самописца. Точное соответствие частоты вращения диска и частоты входного синусоидального сигнала создаёт на втором приборе такую красивую окружность!

Что там за точка M?

Вы заметили на круговом самописце «блуждающую» точку ? Важная, скажу я вам, точка! Это - средняя всех точек, образующих кривую на круговом самописце. Координаты точки получаются как среднее значение -координат и среднее значение -координат всех точек получаемой кривой линии (на примере выше — окружности):

где

- количество точек, образующих кривую линию;
, — и координаты -ой точки кривой .

В физике такая точка называется центром масс механической системы. Наша механическая система состоит из точек единичной массы, образующих в третьем опыте линию окружности (могут быть и другие формы этой кривой). Сами эти точки явным образом изображены на графиках в виде... точек :). Если самописец несколько раз прорисовывает линию (и точки) по одному и тому же пути, мы учитываем все эти точки в наших расчётах.

Более правильное названиедля M — комплексная амплитуда гармоники. В англоязычной литературе встречается слово phasor от phase vector. Соответственно,  — координата M называется вещественной частью комплексной амплитуды, а  — координата называется мнимой частью комплексной амплитуды.

В этой статье комплексные числа не используются, а вместо них применяется геометрическая интерпретация и обычная тригонометрия.

В правом нижнем углу видео отображаются координаты точки M в виде Mx и My, а Md - это длина вектора OM (она же — расстояние от M до начала координат, совпадающего с осью вращения).

Опыт четвёртый: замедляемся...

Давайте замедлим вращение кругового самописца, скажем, в два раза. Пусть время его полного оборота будет 24 с, а сам входной синусоидальный сигнал оставим без изменений:

Кроме того, что на круговом самописце теперь не окружность, а четырёхлепестковый цветочек, замечаем, что по мере накопления данных точка M движется к началу координат (совпадает с осью вращения). Отметим для себя этот момент, а чуть позже вернёмся к нему.

Попробуем замедлить круговой самописец до 1/3 от начальной скорости:

Картинка поменялась, но точка M также, как и в предыдущем примере, предпочитает мигрировать в начало координат.

Опыт пятый: а если крутануть побыстрее?

Разгоним наш круговой самописец до 6 секунд за оборот:

И снова точка M стремится к нулевым координатам!

С чем связано наблюдаемое поведение точки M?

Давайте подумаем, почему при совпадении периодов гармонического сигнала и вращения круга самописца точка M не пытается «прилипнуть» к нулевым координатам, а при несовпадении — прячется в 0.

Итак, у кругового самописца перо перемещается под действием входного сигнала. В нашем случае это синусоида.

Если период колебаний сигнала совпадает с периодом обращения круга, то максимумы и минимумы в положении пера самописца будут приходиться на одни и те же участки (сектора) круга. Например, если максимум перо нарисует при повороте круга на 30 градусов относительно его начального положения, то и при последующих прокрутках круга через этот угол перо будет рисовать максимум. Точно так же, в противоположном секторе, будет повторяться минимум.

Если период колебаний сигнала НЕ совпадает с периодом обращения круга, то максимумы и минимумы будут постоянно сдвигаться в разные положения относительно оси вращения. И в итоге мы получим какую-то достаточно симметричную относительно точки начала координат (0, 0) кривую, у которой M будет приближаться к (0, 0).

Какую информацию можно извлечь из координат точки M?

Интересную. Для сигнала, период которого равен периоду вращения кругового самописца:

• Длина вектора OM равна половине амплитуды или четверти размаха сигнала.

• Угол между осью X на круговом самописце и вектором OM равен фазе сигнала, который в этом случае рассматривается как косинусоидальный. Если мы хотим интерпретировать сигнал как синусоидальный, то следует рассматривать угол между осью Y и OM.

Эта особенность связана с тем, что синус и косинус сдвинуты относительно друга на 90 градусов и иной разницы между ними нет

,

Амплитуда, размах, фаза

Длину вектора OM (на видео длина обозначена Md в правом нижнем углу) можно вычислить по теореме Пифагора, зная координаты точки M:

где

- -координата точки M,
- -координата точки M.

Угол между осью X и вектором OM:

,

Угол между осью Y и вектором OM:

,

Опыт пятый: сдвиг по фазе

Что-то мне надоело следить глазами за крутящимся кругом. Пусть круг будет неподвижен, а вместо этого в обратную сторону вращается перо! В итоге будет такая же картинка, только в более удобном для наблюдения виде.

Рассмотрим сдвинутый по фазе на радиан (30 градусов) синус:

Что мы можем понять из видео? Длина вектора OM стремится к половине амплитуды, то есть, к 5. Угол между вектором OM и осью X стремиться к или 60 градусов — это видно по делениям секторов на круге. Каждый сектор равен 30 градусов. Но, поскольку мы рассматриваем функцию синус, то нас интересует угол между осью Y и вектором OM. А там как раз нужные нам — визуально это ровно 1 сектор на круговом самописце.

Опыт шестой: смешаем краски и разделим их обратно

Попробуем разобрать на составляющие гармоники сложный сигнал, состоящий из трёх компонентов.

Гармоники сигнала — это дополнительные синусоидальные колебания, частоты которых являются кратными основной (самой низкой) частоте сигнала.

Основная гармоника нашего сигнала имеет циклическую частоту радиан/с. Смешаем вместе сигналы 1, 3, 6 гармоник с разными амплитудами и фазами:

Для того, чтобы проанализировать гармоники с 1 по 6 мы должны последовательно посчитать координаты точки M при вращении диска со скоростями, соответствующими гармоникам. Для первой гармоники скорость вращения диска должна быть радиан/с, для второй , для третьей и так далее радиан/с.

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

Первая гармоника:

Видим Md приближается к 1 (половина амплитуды), угол между OM и Y (для синуса) соответствует ожидаемым 60 градусам. Первую гармонику и её параметры мы получить сумели.

Вторая гармоника:

Видно, что M в районе нуля. Делаем вывод, что на второй гармонике тишина.

Третья гармоника:

Md показывает значение, близкое к 2 (половина амплитуды), угол между OM и осью X ожидаемо нулевой, т. к. соответствующий третьей гармонике косинус без начальной фазы.

Четвёртая и пятая гармоники:

Обе этих гармоники показывают отсутствие соответствующих составляющих в исходном сигнале.

Шестая гармоника:

Угол между OM и X около 50 градусов «вниз» (точка M под осью X), что приблизительно соответствует указанной фазе 45 градусов для шестой гармоники, Длина вектора OM около 1.5, что близко к половине амплитуды шестой гармоники.

С увеличением номера (частоты) гармоники точность определения её фазы ухудшается, т.к. в исходном сигнале на короткие периоды высоких гармоник приходится меньше данных (отсчётов), чем на более длительные периоды низкочастотных гармоник.

Это можно заметить визуально на виртуальном круговом самописце. Для этого на графике отображаются «жирные» точки, символизирующие отдельные отсчёты входного сигнала.

Что такое отсчёт?

Отсчёт — численное значение напряжения (или другой измеряемой величины) сигнала в определённый момент времени.

В зарубежной литературе используется термин sample.

На рис. 7. зелёной линией (кривая с точками) изображён график некоего непрерывного сигнала . Дискретные выборки обозначены синими вертикальными линиями. Вот эти выборки с определёнными измеренными значениями сигнала и есть отсчёты или сэмплы. Буквой обозначен период дискретизации или, что тоже самое, период сэмплирования. Естественно, есть и частота дискретизации или частота семплирования, показывающая количество измерений в секунду: .

Опыт седьмой: Приподнятый синус и нулевая гармоника

Рассмотрим функцию:

Что мы можем отметить в этом видео? Хоть это и привычная нам синусоида, но круговой самописец уже не показывает её в виде окружности: «приподнятость» повлияла на форму получаемой кривой.

Длина вектора OM, обозначенная в правом нижнем углу видео как Md, постепенно стремиться к 3, что соответствует половине амплитуды (множитель 6 перед sin). Это ожидаемо. Как и нулевая начальная фаза синусоиды. Всё как если бы мы не прибавляли константу 4. Можно ли на круговом самописце вычислить эту константу?

Можно. Для этого нам нужно построить... нулевую гармонику. Остановить вращение кругового самописца:

Теперь длина OM приближается к 4, что и равно среднему значению функции

Длина OM для нулевой гармоники равна среднему значению всей функции.

Каким образом отсчёты исходного сигнала переносятся на круговой график?

Ранее я упоминал, что мы имеем дело с отдельными отсчётами сигнала. Отсчёт — это мгновенное значение сигнала. Оно получается при измерении исходного аналогового сигнала при помощи специального устройства — аналого‑цифрового преобразователя (АЦП). АЦП выполняет такие измерения с определёнными интервалами времени и мы получаем числа, соответствующие аналоговому сигналу.

В вычислениях используются массивы определённого размера с оцифрованным в виде отдельных отсчётов сигналом. Например, в массив, содержащий 4410 16-битных чисел, можно уместить 0.1 секунду звучания CD Audio одного из двух стерео каналов. Либо 0.05 секунд звучания обоих каналов.

Для отображения этих отсчётов на круговом графике в виде точек кривой нужно знать фазу отображаемого отсчёта и его значение . Со значением всё понятно — это и есть сам отсчёт, а с фазой придётся разобраться.

Фаза — это угол на круговом графике (рис. 8). Угол этот зависит от номера отсчёта в массиве данных из отсчётов и номера гармоники для которой мы строим график.

Если мерить в градусах, то угол равен:

Если в радианах, то:

А почему угол считается именно так?

Давайте рассмотрим случай, когда номер гармоники :

Здесь смысл дроби

в том, чтобы масштабировать позицию отсчёта с номером в массиве (выборке) данных из отсчётов в диапазон значений . Результат этой дроби — множитель, который при умножении на 360 градусов даёт нужный угол . При этом для угол равен 0, а для угол приближается к 360 градусов (он же совпадает с 0 градусов), но не равен ему из-за . Все промежуточные значения равномерно отобразятся в нужные углы .

Коэффициент масштабирует угол для гармоники таким образом, чтобы один период -й гармоники размещался на полном обороте вокруг точки O (рис. 8).

Kоординаты отсчётов сигнала на круговом графике

Координаты, состоящие из угла и радиуса , называются полярными.

Однако, нам нужны привычные декартовы координаты точки . И в этом нет ничего сложного!

Эти формулы следуют из определения синуса и косинуса:

Подставляя формулу для

получим:

где - величина - го отсчёта сигнала, - номер интересующей нас гармоники, - номер отсчёта, для которого мы вычисляем координаты точки .

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

Используя формулу для координат точки :

заменяем и на выражения для и :

где
— вещественная часть комплексной амплитуды  — й гармоники,
 — мнимая часть комплексной амплитуды  — й гармоники,
 — количество отсчётов входного сигнала в выборке для анализа,
 — величина n — го отсчёта,
— номер гармоники от 0 до .

Полученная пара формул и есть дискретное преобразование Фурье (ДПФ).

Расшифровываем телефонные номера из аудиофайла

Наших теоретический знаний уже достаточно для первой практической работы - расшифровки сигналов двухтонального многочастотного набора телефонного номера (Dual-Tone Multi-Frequency, DTMF).

DTMF пришёл на смену импульсному набору как более быстрая и технологичная альтернатива. На постсоветском пространстве DTMF активно внедрялся в сети телефонной связи общего пользования с середины 1990-х годов. В своём городе я впервые услышал его в 1997.

DTMF позволяет закодировать и передать в аналоговую линию связи символы цифр от 0 до 9, дополнительно к ним символы * и # и ещё четыре буквы A. B. C. D. Каждый такой символ передаётся в линию в виде сигналов переменного тока определённых частот. Для формирования одного символа складываются два гармонических сигнала с частотами из таблицы:

А вот на этой страничке можно поиграться (а кому-то и понастальгировать) с Online DTMF генератором.

Для нашего опыта я взял готовую запись очень быстрых телефонных наборов вот отсюда и для удобства работы сконвертировал в WAV формат 8000 сэмплов/с, 8 бит (1 байт) на сэмпл. Готовый WAV-файл в моём git-репозитории.

Код для детектирования телефонных номеров я написал на Python и сейчас мы рассмотрим использованные в нём типичные для подобных задач подходы.

import math as m

# DTMF harmonic resolution.
harmonic_resolution = 35 # 35 Hz per harmonic
# Sample rate.
sample_rate = 8000 # 8000 samples per second
# Harmonic detection magnitude.
harmonic_detection_magnitude = 5
# Calculate the number of samples per chunk.
N = sample_rate // harmonic_resolution # 8000 / 35 = 228 samples per chunk
# Recalculate the harmonic resolution to be the actual value.
harmonic_resolution = sample_rate / N # 8000 / 228 = 35.09 Hz
# Calculate the number of samples per shift.
shift_interval = N // 3 # 228 / 3 = 76 samples per shift
# Silence between whole phone numbers is 0.2 seconds.
silence_interval = 0.2
# Calculate the number of iterations for the silence interval.
silence_number = round(silence_interval / (shift_interval / sample_rate))

# DTMF frequencies. See https://en.wikipedia.org/wiki/Dual-tone_multi-frequency_signaling
dtmf_frequencies = [697, 770, 852, 941, 1209, 1336, 1477, 1633]
# DTMF digits to frequency mapping.
freq_to_digit = {
     (697, 1209): '1', (697, 1336): '2', (697, 1477): '3', (697, 1633): 'A',
     (770, 1209): '4', (770, 1336): '5', (770, 1477): '6', (770, 1633): 'B',
     (852, 1209): '7', (852, 1336): '8', (852, 1477): '9', (852, 1633): 'C',
     (941, 1209): '*', (941, 1336): '0', (941, 1477): '#', (941, 1633): 'D'
}

# Get the harmonic of a frequency.
# frequency is an integer
# Returns the harmonic.
def get_harmonic(frequency: int) -> int:
    return round(frequency / harmonic_resolution)

# DTMF harmonics
dtmf_harmonics = [get_harmonic(f) for f in dtmf_frequencies]

# DTMF digits to harmonics mapping.
harmonic_to_digit = {(get_harmonic(d[0]), get_harmonic(d[1])): h for d, h in freq_to_digit.items()}

# Get the digit from the harmonics.
# harmonic1 and harmonic2 are integers
# Returns the digit or None if the harmonics are not valid DTMF harmonics.
def get_digit(harmonic1: int, harmonic2: int) -> str:
    h = (min(harmonic1, harmonic2), max(harmonic1, harmonic2))
    return harmonic_to_digit.get(h) # get the digit from the mapping


# Get the frequency of the harmonic.
# harmonic is an integer
# Returns the frequency.
def get_frequency(k):
    return k * harmonic_resolution


# Get the magnitude of the harmonic.
# phase_vector is a dictionary with 're' and 'im' keys
# Returns the magnitude.
def get_magnitude(phase_vector: dict) -> float:
    # Magnitude of the harmonic is the square root of the sum of the squares of the real and imaginary parts.
    re = phase_vector['re']
    im = phase_vector['im']
    return m.sqrt(re * re + im * im)


# Calculate the DFT for each harmonic.
# data is a list of samples
# harmonics is a list of harmonics
# Returns a list of phase vectors as {re: float, im: float} for each harmonic.
def calculate_dft(data: list[int], harmonics: list[int]) -> list[dict]:
    dft = [None] * len(harmonics) # prepare the list for {re: float, im: float} for each harmonic
    for i in range(len(harmonics)): # for each harmonic, i is the index
        k = harmonics[i] # k is the harmonic number
        sum_x = 0.0 # sum of the real parts
        sum_y = 0.0 # sum of the imaginary parts
        for n in range(N): # for each sample
            r = data[n] # get the sample
            x = r * m.cos(2 * m.pi * k * n / N) # real part of the sample
            y = r * m.sin(2 * m.pi * k * n / N) # imaginary part of the sample
            sum_x += x # accumulate the real part
            sum_y += y # accumulate the imaginary part
        sum_x /= N # divide the accumulated real part by N to get the average
        sum_y /= N # divide the accumulated imaginary part by N to get the average
        dft[i] = {'re': sum_x, 'im': sum_y} # store the average real and imaginary parts
    return dft


samples = [] # List of samples loaded from the WAV file.
found_dtmf = False # Found a DTMF tone in current chunk.
silence_counter = 0 # Number of iterations with no DTMF tones.

# The WAV file is converted from https://en.wikipedia.org/wiki/File:DTMF_dialing.ogg#file
# It contains 8 phone numbers of 10 DTMF digits each:
# 0696675356, 4646415180, 2336731416, 3608338160, 4400826146, 6253689638, 8482138178, 5073643399
with open("dtmf.wav", "rb") as f:
    # Skip the WAV header, we use the data section which
    # starts at byte 44 and we know it is 8-bit PCM data.
    f.seek(44)
    # Read one chunk and then slide the window by one sample.
    while (byte := f.read(1)): # Read one byte at a time.
        samples.append(byte[0]) # Convert byte to int and append to samples.
        if len(samples) == N:
            # Compute the DFT for each harmonic.
            dft = calculate_dft(samples, dtmf_harmonics)
            magnitudes_per_harmonic = {} # Dictionary of magnitudes per harmonic.
            for i in range(len(dft)): # For each harmonic in the DFT.
                phase_vector = dft[i]
                magnitude = get_magnitude(phase_vector)
                # Store the magnitude of the harmonic.
                magnitudes_per_harmonic[dtmf_harmonics[i]] = magnitude
            # Sort the harmonics by magnitude in descending order and create a list
            # of dictionaries with 'k' and 'm' keys for the harmonic number and magnitude.
            # The first two harmonics are the most significant.
            sorted_harmonics = [
                    {'k': item[0], 'm': item[1]}
                    for item in sorted(magnitudes_per_harmonic.items(), key=lambda item: item[1], reverse=True)[:2]
                ]
            # If the first two harmonics are above the threshold, we probably detected a DTMF tone.
            if len(sorted_harmonics) == 2 and \
                sorted_harmonics[0]['m'] > harmonic_detection_magnitude and \
                sorted_harmonics[1]['m'] > harmonic_detection_magnitude:
                # Try to get the digit from the harmonics.
                digit = get_digit(sorted_harmonics[0]['k'], sorted_harmonics[1]['k'])
                if digit:
                    silence_counter = 0
                    if not found_dtmf:
                        print(digit, end="")
                        found_dtmf = True
                else: # No digit found.
                    print("?", end="")
            else: # The chunk has no DTMF tones.
                found_dtmf = False
                silence_counter += 1

            if silence_counter == silence_number:
                # Print a newline because we have reached the silence interval.
                print("")

            # Shift the window forward by shift_interval samples.
            samples = samples[shift_interval:]

Первым делом нужно определиться с требованиями

Исходные параметры оцифрованного сигнала таковы:

• Частота дискретизации (она же частота сэмплирования): 8000 выборок/с

• Размер отсчёта (сэмпла): 8 бит (1 байт), число без знака.

Обратим внимание на табличку частот DTMF на рис. 10. Ближайшие друг к другу частоты, которые нужно различать в сигнале, 697 и 770 Гц. Между ними 73 Гц. Это минимальное требуемое разрешение по частоте. При таком разрешении эти две гармоники будут соседними и, вполне возможно, что ДПФ покажет некоторую заметную амплитуду на обоих гармониках, даже если на входе сигнал будет содержать только одну из них. Это называется растеканием спектра.

Растекание спектра (Spectral leakage) на соседние гармоники и происходит, когда на длине N массива данных сигнала не помещается целое число периодов исследуемой гармоники.

Для большей надёжности обнаружения нужных нам сигналов выберем вдвое более высокое разрешение по частоте, а именно 35 Гц.

Это разрешение совпадает с первой гармоникой. Частота второй гармоники будет 70 Гц, третьей — 105 Гц и так далее.

Полный период первой гармоники должен соответствовать размеру выборки отсчётов :

где
— частота дискретизации, у нас 8000 сэмплов/с,
 — частота первой гармоники, у нас 35 Гц

Поскольку — это количество отсчётов (сэмплов) во входном массиве (выборке) для ДФТ, то округлим это значение до 228 и пересчитаем частоту первой гармоники :

Получили уточнённую частоту первой гармоники 35.1 Гц.

Всё вышеописанные соображения реализованы в начале программы в нескольких строчках Python-кода:

# DTMF harmonic resolution.
harmonic_resolution = 35 # 35 Hz per harmonic

# Sample rate.
sample_rate = 8000 # 8000 samples per second

# Calculate the number of samples per chunk.
N = sample_rate // harmonic_resolution # 8000 / 35 = 228 samples per chunk

# Recalculate the harmonic resolution to be the actual value.
harmonic_resolution = sample_rate / N # 8000 / 228 = 35.09 Hz

Соответствие частот DTMF номерам наших гармоник

Легко заметить, что частоты DTMF чаще всего не точно соответствуют частоте гармоники. Некоторые частоты попадают в середину между гармониками. Например, 1209 Гц посередине между гармониками 34 и 35. В этом случае ДПФ покажет некоторую амплитуду на обеих гармониках. Это - вторая причина, почему я выбрал требуемое разрешение 35, а не 73 Гц.

С частотами и гармониками разобрались. Теперь подумаем об амплитудах гармоник.

Существуют разные подходы к тому, как определять наличие полезного сигнала на фоне шумов. Мы не будем залезать в дебри. Сэмплы нашего сигнала принимают значения от 0 до 255 со средним значением 128. Вычисление первой и последующих гармоник отбрасывают среднее и остаётся только амплитуда (а точнее, её половина) переменной части синусоиды. Давайте просто «с потолка» примем, что если амплитуда гармоники превышает 10, то она несёт в себе полезный сигнал, а если меньше, то там шумы (или мы «проезжаем» через границу между тональным сигналом и паузой, но об этом позже):

# Harmonic detection magnitude.
harmonic_detection_magnitude = 5

harmonic_detection_magnitude - это половина амплитуды.

Вот и все выкладки относительно исходных параметров вычислений.

Вкусные детальки

Помните, входной блок данных содержит 228 отсчётов? Конечно, чем их больше, тем лучше точность, но...

Подумайте вот о чём: между отсчётами в реальном мире проходит 1/8000 секунды. 228 отсчётов - это 28.5 миллисекунд времени. Чем больше отсчётов нам нужно для вычислений, тем больше времени надо их накапливать во входном буфере ДПФ и тем позже, даже если сам ДПФ работает мгновенно, мы получим результат. В нашем случае все исходные данные кто-то записал в файл и задержка ни на что не влияет, но если вы имеете дело с сигналом в реальном времени, например, вы сделали эквалайзер и хотите управлять звуком, поступающим с музыкальных инструментов, задержка на накопление буфера для ДПФ может повлиять на качество и даже практическую целесообразность вашего решения.

Для детектирования ДПФ важно, чтобы длительность DTMF сигнала не была меньше, чем этот самый входной буфер. По стандарту ETSI ES 201 235-3 для уверенного распознавания длительность DTMF сигнала не должна быть менее 40 мс. То есть, мы с нашими 28.5 мс должны успешно его декодировать.

Входной поток сэмплов наша программа считывает из wav файла. Сначала заполняется входной буфер, а потом, по мере вычисления и анализа гармоник, самые «старые» сэмплы из начала буфера удаляются, а в конец дочитывается столько же новых сэмплов из файла и снова запускается ДПФ. Зачем так сделано? почему бы не читать последовательно блоки данных размером по 228 байт и сразу обрабатывать? Дело в том, что искомый сигнал не обязательно начнётся с начала блока. Может и с середины, например. Поэтому я решил сдвигать это 228-байтовое окошко на 1/3 его ширины, то есть, на 76 байт. Таким образом, где бы сигнал DTMF ни начался, мы уверенно его обнаружим:

# Calculate the number of samples per shift.
shift_interval = N // 3 # 228 / 3 = 76 samples per shift

После каждого DTMF символа следует такая же по длительности пауза. Зачем я об этом пишу? Дело в том, что по мере сдвига «окна» входного буфера на одну треть его размера, мы вполне можем детектировать один и тот же DTMF символ несколько раз за несколько последовательных сдвигов. Поэтому, решение о том, что это не один символ «4», а два подряд идущих символа «44» мы можем принять только после обнаружения такой же по длительности тишины на DTMF гармониках между символами. Этот момент учитывается в кусочке:

                    if not found_dtmf:
                        print(digit, end="")
                        found_dtmf = True

Выводим на консоль обнаруженный символ, если ранее флаг обнаружения found_dtmf был сброшен. Затем устанавливаем флаг обнаружения.

Сбрасывается флаг обнаружения только при детектировании тишины на DTMF гармониках:

            else: # The chunk has no DTMF tones.
                found_dtmf = False

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

Запускаем

Программа читает dtmf.wav и выдаёт построчно обнаруженные там DTMF-кодированные номера телефонов.

Заключение

Я надеюсь, что вы найдёте эту статью полезной. Многие свойства дискретного преобразования Фурье я не упомянул и не ставил перед собой такой задачи. Работа над симулятором самописцев и примером заняла в сумме около двух недель по вечерам и ночам после основной работы. Я попытался максимально доходчиво, в картинках, показать как работает ДПФ и где его можно применить. В угоду наглядности отброшены все строгие математические выкладки и выводы.

Буду благодарен конструктивным комментариям и замечаниям. Код здесь.

А. Галилов