Как я решила уравнение жвачки

Однажды на просторах интернета мне попался ролик о связи физики и духовных практик. Там прозвучала короткая фраза:

Ну сейчас пока мы не можем, жуя резинку, решить уравнение.

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

Хорошо, решить уравнение, жуя жвачку, мы пока не можем. А что, если попробовать решить уравнение самой жвачки?

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

Так появился маленький проект на Python, в котором я взяла запись жевания жвачки и попыталась получить из неё уравнение.

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

Что значит «превратить звук в уравнение»

В цифровой записи звук уже представлен числами. Через одинаковые промежутки времени программа сохраняет значение амплитуды звуковой волны.

Упрощённо это выглядит так:

время: 0.0000 0.0001 0.0002 0.0003 ...

амплитуда: 0.012 -0.008 0.031 0.027 ...

Получается набор точек (t_i, y_i), где:

t_i — время;
y_i — амплитуда звука в этот момент.

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

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

В общем виде моя формула выглядит так:

$x(t) \approx \mu + \sum_{k=1}^{N} A_k \cos(2\pi f_k t + \varphi_k)$

Здесь:

x(t) — амплитуда звука в момент времени $t$;
$\mu$ — среднее значение сигнала;
A_k — амплитуда отдельной гармоники;
f_k — её частота;
$\varphi_k$ — фаза;
N — количество гармоник, которые мы решили оставить.

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

Что было на входе

Я нашла в интернете короткую запись жевания жвачки. Это MP3-файл:

длительность — примерно 6,06 секунды;
один аудиоканал;
исходная частота дискретизации — 44 100 Гц.

44 100 Гц означает, что каждую секунду в записи хранится 44 100 значений амплитуды. За шесть секунд получается больше 267 тысяч отсчётов.

Для эксперимента такая точность была избыточной, поэтому перед анализом я уменьшила частоту дискретизации до 8000 Гц. Этого достаточно, чтобы увидеть структуру сигнала, а вычисления и итоговые файлы становятся заметно компактнее.

Сверху на графике видна форма записи во времени, снизу — её частотный спектр. Отдельные сильные всплески на верхнем графике — как раз те самые щелчки и резкие звуки, из-за которых жевание сложно описать небольшой красивой формулой.

Шаг 1. Декодируем MP3

Сам Python не обязан уметь читать все аудиоформаты, поэтому для декодирования я использовала ffmpeg. На выходе получался монофонический поток 32-битных чисел с плавающей точкой:

command = [ “ffmpeg”, “-v”, “error”, “-i”, str(path), “-ac”, “1”, “-ar”, str(sample_rate), “-f”, “f32le”, “-acodec”, “pcm_f32le”, “-”, ]

result = subprocess.run(command, check=True, stdout=subprocess.PIPE) samples = np.frombuffer(result.stdout, dtype=np.float32)

Затем я нормализовала сигнал, чтобы все значения находились примерно в диапазоне от -1 до 1:

peak = np.max(np.abs(samples))

if peak > 0: samples = samples / peak

После этого звук действительно превращается в обычный массив чисел NumPy.

Шаг 2. Раскладываем звук на частоты

Чтобы найти гармоники, я использовала быстрое преобразование Фурье, или FFT:

centered = samples - np.mean(samples) spectrum = np.fft.rfft(centered) frequencies = np.fft.rfftfreq( centered.size, d=1.0 / sample_rate, ) amplitudes = (2.0 / centered.size) * np.abs(spectrum) phases = np.angle(spectrum)

После этого для каждой найденной частоты у нас есть три нужных числа:

частота + амплитуда + фаза

Именно из них затем собираются слагаемые вида:

$A_k\cos(2\pi f_k t+\varphi_k)$

Шаг 3. Выбираем главные компоненты

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

Поэтому я оставила 32 заметные спектральные компоненты. При этом соседние частоты должны отличаться хотя бы на 25 Гц. Без этого почти все места в формуле занимали бы очень близкие частоты вокруг одного сильного пика.

Код отбора выглядит примерно так:

strongest_first = usable[np.argsort(amplitudes[usable])[::-1]]

selected = []

for index in strongest_first:

frequency = frequencies[index]

if all(abs(frequency - frequencies[chosen]) >= 25 for chosen in selected ): selected.append(index)

if len(selected) == 32: break

Это уже компромисс. Чем больше компонентов оставить, тем ближе восстановленный сигнал будет к исходному. Но тем длиннее и страшнее станет формула.

Какое уравнение получилось

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

$x(t) \approx \mu + \sum_{k=1}^{32} A_k\cos(2\pi f_k t+\varphi_k), \qquad 0 \leq t \leq 6{,}0605$

где:

$\mu = 3{,}7733777 \cdot 10^{-5}$

Все 32 набора коэффициентов я сохранила в CSV. Но для статьи гораздо приятнее выглядит сокращённая версия из восьми самых сильных компонентов:

$\begin{aligned} x_8(t) \approx 10^{-3}\big(&0{,}038 \\ &+1{,}917\cos(2\pi\cdot885{,}57t+0{,}156)\\ &+1{,}559\cos(2\pi\cdot56{,}76t+2{,}718)\\ &+1{,}537\cos(2\pi\cdot912{,}47t-0{,}461)\\ &+1{,}504\cos(2\pi\cdot859{,}34t-0{,}567)\\ &+1{,}133\cos(2\pi\cdot827{,}82t+1{,}413)\\ &+1{,}007\cos(2\pi\cdot30{,}69t-1{,}690)\\ &+0{,}990\cos(2\pi\cdot797{,}95t-2{,}434)\\ &+0{,}696\cos(2\pi\cdot937{,}55t-2{,}979)\big) \end{aligned}$

Вот это я и называю уравнением жвачки.

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

Можно ли услышать результат

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

t = np.arange(samples_count) / sample_rate

y = np.full(samples_count, bias)

for term in terms:

y += term.amplitude np.cos( 2.0 np.pi term.frequency_hz t + term.phase_rad )

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

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

А можно описать звук точнее?

Да. Самый прямой способ — оставить все отсчёты и задать между соседними точками линейную функцию:

$x(t)=y_i+ \frac{y_{i+1}-y_i}{t_{i+1}-t_i}(t-t_i), \qquad t_i\leq t\leq t_{i+1}$

Это кусочно-линейная функция. Для моей записи при частоте 8000 Гц получилось около 48 тысяч точек.

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

Поэтому у проекта получилось два ответа:

Короткое и красивое уравнение Фурье — хорошо показывает идею и основные частоты.
Кусочно-линейная функция по всем точкам — лучше описывает конкретную запись, но выглядит как большой CSV-файл.

Что в итоге

Конечно, я не нашла универсальное уравнение жевания жвачки. Такого одного уравнения в рамках этого эксперимента и быть не могло.

Я получила математическую модель конкретного звукового фрагмента. Причём в двух вариантах: компактном, через сумму гармоник, и более точном, через набор точек.

Но мне в этой истории интереснее даже не результат, а сама цепочка:

случайная фраза -> странный вопрос -> звук -> массив чисел -> преобразование Фурье -> коэффициенты -> уравнение -> звук, восстановленный из уравнения

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

Решить уравнение, жуя жвачку, у меня пока не получилось. Зато получилось решить уравнение жвачки.

ссылка на оригинал статьи https://habr.com/ru/articles/1044610/