Документация Engee
Notebook

Спектральный анализ и фильтрация аудиофайла

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

Подключение библиотек и визуализация аудио-данных

Нам понадобятся основные библиотеки для цифровой обработки сигналов и работы с WAV-файлами:

In [ ]: 
#Pkg.add("DSP")
#Pkg.add("WAV")
In [ ]: 
using DSP
using WAV
using Base64

Мы импортируем аудио-файл, содержащий смесь голосового сообщения (обратный отсчёт), и наложенный на него шум толпы, предварительно пропущенный через фильтр верхних частот. Подобный синтезированный сигнал позволит более наглядно разделить частотные диапазоны полезного сигнала и шума.

In [ ]: 
filename = "$(@__DIR__)/countdown_noisy.wav"
audio_data, fs = wavread(filename);
nsamples = length(audio_data);
dt = 1/fs;
t = 0:dt:((nsamples - 1) .* dt);
original = audio_data[1:nsamples];

Здесь мы объявляем пользовательскую функцию, позволяющую прослушивать WAV-файлы и использовать интерактивный инструмент проигрывателя:

In [ ]: 
function audioplayer(s, fs);
    buf = IOBuffer();
    wavwrite(s, buf; Fs=fs);
    data = base64encode(unsafe_string(pointer(buf.data), buf.size));
    markup = """<audio controls="controls" {autoplay}>
                <source src="data:audio/wav;base64,$data" type="audio/wav" />
                Your browser does not support the audio element.
                </audio>"""
    display("text/html", markup);
end

Прослушаем исходный аудио-файл (полезный сигнал на фоне высокочастотного шума):

In [ ]: 
audioplayer(original, fs)

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

In [ ]: 
plot(t[1:10:end], original[1:10:end], 
    xlabel = "Время, с", 
    ylabel = "Амплитуда", 
    title = "Исходный сигнал во времени")
Out[0]:

Спектральный анализ

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

In [ ]: 
endindx = 2^16;
sig_crop = original[1:endindx];
p1 = DSP.welch_pgram(sig_crop, onesided=true, fs=fs);
plot(freq(p1), pow2db.(power(p1)), 
    xguide = "Частота, Гц", 
    yguide = "Мощность/частота, дБ/Гц",
    title = "Спектр исходного сигнала")
Out[0]:

Шум нестационарен во времени, поэтому полезно оценить и спектрограмму сигнала - то есть зависимость спектра от времени. Для этого воспользуемся функцией spectrogram:

In [ ]: 
spgr = DSP.spectrogram(original, 511, 256; fs = fs);
heatmap(spgr.time, spgr.freq, pow2db.(spgr.power), 
        ylim=(0,10000), 
        clim=(-140,-40), 
        xlabel = "Время, сек", 
        ylabel = "Частота, Гц",
        title = "Спектрограмма исходного сигнала")
Out[0]:

Как видно из спектрограммы, на всём протяжении аудиофайла шум в спектре занимает полосу выше 3 кГц, а это значит, что речевой сигнал можно выделить стандартным фильтром нижних частот (ФНЧ).

Разработка цифрового фильтра

Разработаем цифровой фильтр с подходящими характеристиками. Для создания фильтра нижних частот мы используем функцию digitalfilter, в которой мы указываем тип отклика фильтра lowpass с частотой пропускания 3000 Гц, а также объявляем метод его синтеза Elliptic: это соответствует эллиптическому фильтру с бесконечной импульсной характеристикой (БИХ) и требуемыми характеристиками (порядок фильтра - 10, пульсации в полосе пропускания - не более 0.1 дБ, подавление в полосе заграждения - не ниже 60 дБ):

In [ ]: 
myfilt = digitalfilter(Lowpass(3000; fs), Elliptic(10, 0.1, 60));

Рассчитаем и отобразим характеристики фильтра. В первую очередь нас интересует амплитудно-частотная характеристика (АЧХ), позволяющая понять, какие частоты спектра исходного сигнала наш фильтр будет пропускать, а какие - подавлять. Для этого воспользуемся функцией freqresp:

In [ ]: 
H, w = freqresp(myfilt);
freq_vec = fs*w/(2*pi);
plot(freq_vec, pow2db.(abs.(H))*2, linewidth=3,
     xlabel = "Частота, Гц",
     ylabel = "Амплитуда, дБ",
     title = "АЧХ фильтра нижних частот")
Out[0]:

Теперь отобразим фазо-частотную характеристику (ФЧХ) фильтра функцией phaseresp:

In [ ]: 
phi, w = phaseresp(myfilt);
plot(freq_vec, rad2deg.(phi/2), linewidth = 3, 
    xlabel = "Частота, Гц",
    ylabel = "Фаза, град",
    title = "ФЧХ фильтра нижних частот")
Out[0]:

Наконец, рассчитаем и отобразим импульсную характеристику (ИХ) фильтра (точнее, первых её 100 отсчётов, так как наш фильтр имеет бесконечную ИХ) функцией impresp:

In [ ]: 
myimpresp = impresp(myfilt,100);
plot(myimpresp, line=:stem, marker=:circle, 
    xlabel = "Отсчёты ИХ", 
    ylabel = "Значения ИХ",
    title = "Импульсная характеристика ФНЧ")
Out[0]:

Применим разработанный фильтр к исходному аудио-сигналу функцией filt и отобразим результат обработки в частотной области:

In [ ]: 
filtered = filt(myfilt, original);
outsig = filtered[1:endindx];
p2 = DSP.welch_pgram(outsig, onesided=true, fs=fs);
plot(freq(p1), [pow2db.(power(p1)), pow2db.(power(p2))], 
    xlabel = "Частота, Гц", 
    ylabel = "Мощность/частота, дБ/Гц",
    title = "Спектр обработанного сигнала")
Out[0]:

На спектре можно убедиться, что частоты выше 3 кГц были успешно подавлены.

Также отобразим результат фильтрации во временной области:

In [ ]: 
plot(t[1:10:end], [original[1:10:end], filtered[1:10:end]], 
    xlabel = "Время, с", 
    ylabel = "Амплитуда", 
    title = "Обработанный сигнал во времени")
Out[0]:

Из графика видно, что между произносимыми словами нет шумового заполнения. Убедимся в успешности обработки, прослушав отфильтрованный аудио-сигнал:

In [ ]: 
audioplayer(filtered, fs)

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

Заключение

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