Чем меньше битрейт потока данных, тем дальше можно его передать. Если получится сжать голос до 1-2 кБ/с, можно будет общаться на расстоянии до 25км (с использованием трансивера CC1125). Давайте попробуем разработать голосовой кодек самостоятельно, без использования готовых решений вроде codec2.

Как устроен голос? В нём можно выделить несколько элементов:

• гласные звуки
• согласные — шумовые — взрывные: п, к, т…
• согласные — шумовые — фрикативные/свистящие: с, ш, щ, ф, в…
• согласные — сонорные — дрожащие: р
• согласные — сонорные — носовые: м, н

Гласные звуки можно рассматривать как музыкальные — это чистый тон + его обертоны различной амплитуды. Согласные звуки — это шумовые звуки с разным тембром и сонорные (тон + шум).

Гласные и согласные здорово различаются по своей информативности, сравните:

гласные: .и..ио.ы .и..ио.ы .и..ио.ы а.ы. .о.

согласные: м.лл..н. м.лл..н. м.лл..н. .л.х р.з

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

Давайте начнём с гласных, раз уж их так легко описать. Я разрабатывал алгоритм на Питоне, просто потому что в нём это реально быстро и просто.

Чтение звукового файла

Я записал звук «а» в собственном исполнении и в исполнении девушки, и сохранил его в сыром виде с дикретизацией 8кГц в виде 16-битных знаковых целых с помощью Audacity. Сделаем функцию для чтения такого файла:

from struct import *

def read_file(filename):
f = open(filename, "rb")
chunk = f.read(2)
raw_input = []
while chunk:
raw_input.append(unpack('<h', chunk)[0])
chunk = f.read(2)
f.close()
return raw_input

raw_input = read_file("a_cat_m_8k_200.raw")[00]

Я сохраняю только 800 отсчётов, т.е. 100мс звука — на таком отрезке ещё можно принять сигнал стационарным, неизменяющимся. Следующие фреймы потом проанализируем очередными проходами алгоритма. Построим осциллограмму полученного звука:

import matplotlib.pyplot as plt
plt.plot(cepstrum[0:100], 'r')
plt.grid(b=True, which='major', color='b', linestyle='-')
plt.show()

Построение спектра

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

import cmath

def dft(fnList):
pi2 = cmath.pi * 2.0
N = len(fnList)
FmList = []
for m in range(N):
Fm = 0.0
for n in range(N): Fm += fnList[n] * cmath.exp(- 1j * pi2 * m * n / N)
FmList.append(abs(Fm / N))
return FmList

spectrum = dft(raw_input)

plt.plot(spectrum[0:400], 'r')

Красивый спектр, правда? Увеличим интересующий нас фрагмент:

Хорошо виден основной тон и его обертоны, спектр имеет периодическую структуру.

Положение каждого пика определяется соотношением: freq = samples / period. Абсцисса пика равна отношению длины фрейма в семплах (в нашем случае 800) к периоду волны.

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

Чтобы уменьшить этот эффект, необходимо применить к сигналу так называемое «окно». Я умножу сигнал на окно Хэннинга.

raw_input = [amp * 1.8 * (0.5 - 0.5 * math.cos(2*3.1415*i/len(raw_input))) for i, amp in enumerate(raw_input)]

Красная линия — спектр исходного сигнала, зелёная — спектр сигнала с окном Хэннинга:

Стало гораздо лучше: пики выделены сильнее, их полуширина меньше и в спектре меньше шума.

Построение кепстра

Здесь уже видно всё, что нам нужно, но как программно найти пики и их амплитуды? Проходить обычным поиском локального максимума — неблагодарное дело: пики могут быть слегка сдвинуты, дублированы, и так мы найдём только целую часть периода, а он может быть дробным. Как же нам тогда определить этот период? Тем же способом, как мы определяли периоды в исходном звуковом сигнале, преобразованием Фурье. Только нужно сначала логарифмировать наш спектр, чтобы сильнее выделить пики.

И применить к этому преобразование Фурье:

def calc_cepstrum(spectrum):
return dft([math.log(amp) for i, amp in enumerate(spectrum)])
cepstrum = calc_cepstrum(spectrum)
plt.plot(cepstrum[0:100], 'b')

То, что мы построили, называется прикольным словом «кепстр» (Cepstrum) — из названия понятно, что это является чем-то вроде обращения спектра.

Поиск фундаментальной частоты

В нём нам нужен только первый пик, чьё положение будет в точности равно периоду искомого основного тона. Мы найдём этот пик обычным поиском глобального максимума, ограничив область поиска частотами 80..270 Гц, т.к. именно в этом диапазоне лежат основные частоты голосов взрослых людей. В нашем случае это диапазон периодов 30..100 (8000/270..8000/80).

import numpy as np
fund_freq = len(cepstrum) / (np.argmax(cepstrum[30:100]) + 30.0)
print 8000 * fund_freq / len(cepstrum)

> 112.676

В моём случае фундаментальная частота была равна 113Гц, а у моей девушки — ровно 200Гц, что хорошо соотносится с цитатой из Википедии:

Голос типичного взрослого мужчины имеет фундаментальную частоту (нижнюю) от 85 до 155 Гц, типичной взрослой женщины от 165 до 255 Гц.

Запись обертонов

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

harms = []
for i in range(1, 10): harms.append(int(sum([spectrum[int(round(i * fund_freq)) + j] for j in range (-1, +1) ])))
print harms

> [952, 915, 1111, 674, 803, 898, 988, 879, 463]

Даже девятая гармоника имеет амплитуду, сравнимую с амплитудой основной частоты, поэтому нужно записать довольно много обертонов, минимум 8-10. С другой стороны, нам достаточно ограничиться полосой в 4кГц, более высокие гармоники уже не так нужны для передачи голоса.

Упаковка параметров

Для описания одного фрейма с голосом (10мс в нашем случае) требуется:

• фундаментальная частота — 2 байта
• N * амплитуда обертона — N/2 байт

Обычно достаточно не более 10 обертонов, то есть характеристика займёт 2 + 5 = 7 байт. Исходный фрейм весил 1600 байт, то есть степень сжатия составляет ~230 раз. Но нужно понимать, что так мы сжали только гласные звуки, ещё остались согласные — которые, как назло, передают куда больше информации.

Распаковка звука и сравнение

Разжатие звука обратно можно было бы сделать через обратное преобразование Фурье, но сделаем ещё проще — сгенерируем синусоиды нужных частот и амплитуд, и сложим их все друг с другом.

out = []
for i in range(0, 800):
out.append((harms[0] * math.sin(2 * 3.1415926 * i * 1 / (len(cepstrum) / fund_freq)) +
harms[1] * math.sin(2 * 3.1415926 * i * 2 / (len(cepstrum) / fund_freq)) +
harms[2] * math.sin(2 * 3.1415926 * i * 3 / (len(cepstrum) / fund_freq)) +
harms[3] * math.sin(2 * 3.1415926 * i * 4 / (len(cepstrum) / fund_freq)) +
harms[4] * math.sin(2 * 3.1415926 * i * 5 / (len(cepstrum) / fund_freq)) +
harms[5] * math.sin(2 * 3.1415926 * i * 6 / (len(cepstrum) / fund_freq)) +
harms[6] * math.sin(2 * 3.1415926 * i * 7 / (len(cepstrum) / fund_freq))) * 3)
# print out[0:100]
f = open("a_cat_m_8k_200.out", "wb")
for i in range(0, 800):
f.write(pack('<h', out[i]))
f.close()

Откроем в Audacity исходный файл и восстановленный, и сравним осциллограммы:

И спектры:

Довольно похоже. На слух тоже похоже, причём 10 гармоник звучат гораздо лучше, чем 8.