Александр Иванов – Цифровая обработка сигналов на Python. От инженера к разработчику. (страница 5)

demonstrate_aliasing(signal_freq=25000, sample_rate=44100) # алиасинг!

demonstrate_aliasing(signal_freq=40000, sample_rate=44100) # сильный алиасинг

Разберём, что делает этот код. Функция demonstrate_aliasing принимает частоту сигнала и частоту дискретизации. Она вычисляет частоту Найквиста — половину частоты дискретизации. Если частота сигнала ниже или равна частоте Найквиста, всё в порядке. Если выше — вычисляется частота алиаса по формуле: abs(signal_freq - sample_rate * round(signal_freq / sample_rate)). Эта формула показывает, на какой частоте появится ложный сигнал.

Первый пример: сигнал 5 000 Гц, дискретизация 44 100 Гц. Частота Найквиста — 22 050 Гц. Сигнал ниже, алиасинга нет. Второй пример: сигнал 25 000 Гц при той же дискретизации. Частота сигнала выше частоты Найквиста. По формуле алиас будет на частоте 19 100 Гц. То есть мы пытаемся записать звук частотой 25 000 Гц, а после дискретизации слышим 19 100 Гц — совсем другой тон! Третий пример: сигнал 40 000 Гц. Алиас будет на частоте 4 100 Гц — высокий тон превратился в низкий.

Запустите скрипт и прослушайте файлы. Первый файл звучит как высокий, но чистый тон. Второй — тон заметно ниже, хотя мы генерировали более высокую частоту. Третий — тон стал басовитым, что совершенно не соответствует исходным 40 000 Гц. Это алиасинг в действии.

Борьба с алиасингом: антиалиасинговый фильтр

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

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

python

from scipy.signal import butter, sosfilt

def anti_alias_filter(y, sr, cutoff=None):

"""

Применяет антиалиасинговый фильтр (фильтр нижних частот).

cutoff - частота среза. По умолчанию = частота Найквиста * 0.9

"""

nyquist = sr / 2

if cutoff is None:

cutoff = nyquist * 0.9 # небольшой запас

# Создаём фильтр Баттерворта 8-го порядка

sos = butter(8, cutoff / nyquist, btype='lowpass', output='sos')

y_filtered = sosfilt(sos, y)

return y_filtered

# Демонстрируем работу фильтра

sr_high = 441000

t = np.linspace(0, 1.0, int(sr_high * 1.0), endpoint=False)

# Создаём сигнал с частотами 5000 и 25000 Гц одновременно

y_mixed = np.sin(2 * np.pi * 5000 * t) + 0.5 * np.sin(2 * np.pi * 25000 * t)

# Дискретизируем без фильтра

step = sr_high // 44100

y_bad = y_mixed[::step]

# Применяем антиалиасинговый фильтр, затем дискретизируем

y_filtered = anti_alias_filter(y_mixed, sr_high, cutoff=20000)

y_good = y_filtered[::step]

sf.write('aliasing_bad.wav', y_bad, 44100)

sf.write('aliasing_good.wav', y_good, 44100)

print("Файлы сохранены. Сравните aliasing_bad.wav и aliasing_good.wav")

print("В плохом файле слышен алиас на ~19100 Гц от сигнала 25000 Гц")

print("В хорошем файле сигнал 25000 Гц подавлен фильтром до дискретизации")

Этот код создаёт смесь двух частот: 5 000 Гц (безопасная) и 25 000 Гц (опасная, выше частоты Найквиста для 44 100 Гц). Без фильтра после дискретизации мы услышим алиас от 25 000 Гц на частоте около 19 100 Гц. С фильтром частота 25 000 Гц подавляется до дискретизации, и мы слышим только чистый тон 5 000 Гц. Прослушайте оба файла и убедитесь сами.

Квантование: превращаем высоту волны в число

Дискретизация решает проблему времени: мы знаем, как часто измерять сигнал. Но есть ещё проблема амплитуды. Звуковая волна непрерывна не только во времени, но и по амплитуде. В каждый момент времени звуковое давление может принимать любое значение в некотором диапазоне. Компьютер не может хранить «любое» значение — он может хранить только конечное количество различных чисел. Процесс округления измеренного значения до ближайшего допустимого числа называется квантованием.

Представьте, что вы измеряете рост человека. Вы можете записать «175 сантиметров», но не «175,384927... сантиметров». Вы округляете до целых сантиметров. Это и есть квантование. В цифровом аудио амплитуда округляется до ближайшего значения, которое может быть представлено заданным количеством бит.

Количество бит определяет, сколько различных уровней громкости мы можем различить. При 8 битах у нас есть 2 в 8-й степени = 256 различных уровней. При 16 битах — 65 536 уровней. При 24 битах — более 16 миллионов уровней. Чем больше бит, тем точнее мы можем записать амплитуду и тем меньше шум квантования — разница между истинным значением и его квантованным приближением.

Шум квантования звучит как тихое шипение, добавленное к сигналу. При 16 битах этот шум очень тихий — около -96 децибел относительно максимального сигнала. При 8 битах шум отчётливо слышен. Давайте смоделируем квантование и послушаем разницу.

python

def quantize(y, bits):

"""

Квантует сигнал до указанного количества бит.

bits - количество бит (например, 8 или 16)

"""

levels = 2 ** bits

# Диапазон от -1 до 1 делим на levels уровней

y_normalized = y / np.max(np.abs(y)) # нормализуем к диапазону -1..1

y_quantized = np.round(y_normalized * (levels / 2 - 1)) / (levels / 2 - 1)

return y_quantized

# Загружаем реальный голос

import librosa

y, sr = librosa.load('voice_sample.wav', sr=44100)

# Квантуем с разной разрядностью

for bits in [16, 12, 8, 6, 4]:

y_q = quantize(y, bits)

sf.write(f'voice_{bits}bit.wav', y_q, sr)

# Вычисляем шум квантования

noise = y - y_q[:len(y)]

noise_power = np.mean(noise ** 2)