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
Out[0]:
audioplayer (generic function with 1 method)

让我们来听听原始音频文件（在高频噪音背景下的有用信号）：

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 kHz 以上的频段，这意味着语音信号可以通过标准低通滤波器（LF）分离出来。

数字滤波器设计

让我们设计一个具有合适特性的数字滤波器。要创建一个低通滤波器，我们可以使用函数digitalfilter ，其中指定通带频率为 3000 Hz 的滤波器响应类型lowpass ，并声明其合成方法Elliptic ：这相当于一个具有无限脉冲响应 (IIR) 和所需特性（滤波器阶数 - 10，通带纹波 - 不大于 0.1 dB，栅栏带抑制 - 不小于 60 dB）的椭圆滤波器：

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

让我们来计算并显示滤波器的特性。首先，我们对幅频响应（AFR）很感兴趣，通过它可以了解滤波器将通过原始信号频谱中的哪些频率，抑制哪些频率。为此，让我们使用函数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 来显示滤波器的相频响应（FFR）：

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

最后，使用函数impresp 计算并显示滤波器的脉冲响应（IR）（更准确地说，是前 100 个采样，因为我们的滤波器具有无限的 IR）：

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 kHz 以上的频率已被成功抑制。

我们也来展示一下时域滤波的结果：

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

从图中可以看出，在说话词语之间没有噪音填充。让我们来听听滤波后的音频信号，确保处理成功：

In [ ]: 
audioplayer(filtered, fs)

所应用的滤波器不仅成功去除了信号中的噪音，还抑制了语音信号中的高频成分，因此听起来没有那么明亮的泛音和嘶嘶声。

结论

以选定的音频文件为例，我们了解了 Engee 中典型的数字信号处理工作流程，包括导入、可视化、分析和处理。