音频文件的频谱分析和过滤

在这个例子中，我们考虑一个典型的数字信号处理任务的技术计算工作流程。作为示例的一部分，导入、可视化和收听包含有用信号和高频干扰的混合的音频文件。接下来，在频域中分析信号，选择合适的数字滤波器的规格，执行其合成，并可视化滤波器特性。结果，将开发的滤波器应用于原始信号，并分析和可视化滤波结果。

连接库和可视化音频数据

我们将需要用于数字信号处理和处理WAV文件的基本库。:

#Pkg.add("DSP")
#Pkg.add("WAV")

using DSP
using WAV
using Base64

我们正在导入一个音频文件，其中包含叠加在其上的语音消息（倒计时）和人群噪音的混合，以前通过高通滤波器过滤。这样的合成信号将使得能够更清楚地分离有用信号和噪声的频率范围。

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文件并使用交互式播放器工具。:

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

audioplayer (generic function with 1 method)

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

audioplayer(original, fs)

让我们在时域中显示信号：在图形上，您可以观察到对应于倒计时语音消息的幅度的特征周期性突发（对于时域渲染的速度，细化长期信号是有用的-在我

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

光谱分析

在时域中可视化信号并不能让我们清楚地了解如何将语音信号与噪声分开。因此，有必要对信号进行频谱分析，即评估其频谱（频谱功率密度）。为此，我们将使用该函数 welch_pgram. 值得注意的是，为了节省计算时间，我们只取原始信号的前2^16个样本。

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 = "Спектр исходного сигнала")

噪声在时间上不稳定，因此评估信号的频谱图-即频谱对时间的依赖性是有用的。为此，请使用函数 spectrogram:

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 = "Спектрограмма исходного сигнала")

从频谱图可以看出，在整个音频文件中，频谱中的噪声占据了3kHz以上的频带，这意味着可以用标准低通滤波器（LPF）隔离语音信号。

数字滤波器的开发

我们将开发具有合适特性的数字滤波器。要创建一个低通滤波器，我们使用函数 digitalfilter 其中我们指定过滤器响应的类型 lowpass 传输频率为3000赫兹,我们还宣布了一种合成方法. Elliptic 这对应于具有无限脉冲响应（IIR）和所需特性（滤波器阶数-10，通带中的纹波-不超过0.1db，振臂带中的抑制-不低于60dB）的椭圆滤波器:

myfilt = digitalfilter(Lowpass(3000), Elliptic(10, 0.1, 60); fs);

计算并显示滤波器特性。首先，我们对振幅-频率响应（AFC）感兴趣，这使我们能够了解滤波器将通过原始信号频谱的哪些频率，哪些频率将被抑制。为此，请使用函数 freqresp:

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

现在我们将显示滤波器的相位-频率响应（FCH）与功能 phaseresp:

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

最后，我们将通过函数计算并显示（他们的）滤波器的脉冲响应（更确切地说，它的前100个样本，因为我们的滤波器有无限多个） impresp:

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

让我们将开发的滤波器应用于原始音频信号。 filt 我们将在频域中显示处理的结果。:

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 = "Спектр обработанного сигнала")

在频谱上，可以验证3kHz以上的频率已被成功抑制。

我们还将在时域中显示过滤结果。:

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

从图中可以看出，口语之间没有噪声填充。让我们通过收听滤波后的音频信号来验证处理的成功。:

audioplayer(filtered, fs)

所应用的滤波器成功地从信号中去除了噪声，但也抑制了语音信号的高频分量，这使得它听起来不那么明亮地带有泛音和在聆听时发出嘶嘶声。

结论

使用所选音频文件的示例，我们检查了Engee中数字信号处理的典型工作流程，其中包括导入，可视化，分析和处理。