Engee documentation
Notebook

Spectral analysis and filtering of an audio file

This example discusses a typical workflow of technical computing for digital signal processing tasks. As part of the example, an audio file containing a mixture of a useful signal and high-frequency interference is imported, visualized, and listened to. Next, the signal is analyzed in the frequency domain, the specification of a suitable digital filter is selected, its synthesis is performed, and the filter characteristics are visualized. As a result, the developed filter is applied to the original signal, and the filtering result is analyzed and visualized.

Connecting libraries and visualizing audio data

We will need basic libraries for digital signal processing and working with WAV files.:

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

We are importing an audio file containing a mix of a voice message (countdown) and crowd noise superimposed on it, previously filtered through a high-pass filter. Such a synthesized signal will make it possible to more clearly separate the frequency ranges of the useful signal and noise.

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];

Here we are announcing a custom function that allows you to listen to WAV files and use the interactive player tool.:

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)

Let's listen to the original audio file (a useful signal against a background of high-frequency noise):

In [ ]: 
audioplayer(original, fs)

Let's display the signal in the time domain: on the graph, you can observe characteristic periodic bursts of amplitude corresponding to the countdown voice message (for the speed of rendering in the time domain, it is useful to thin out long-term signals - in our case, we take every 10th countdown).

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

Spectral analysis

Visualizing a signal in the time domain does not give us a clear understanding of how to separate a voice signal from noise. Therefore, it is necessary to perform a spectral analysis of the signal, that is, to evaluate its spectrum (spectral power density). To do this, we will use the function welch_pgram. It is worth noting that to save calculation time, we take only the first 2^16 samples of the original signal.

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]:

Noise is unstable in time, so it is useful to evaluate the spectrogram of the signal - that is, the dependence of the spectrum on time. To do this, use the function 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]:

As can be seen from the spectrogram, throughout the audio file, the noise in the spectrum occupies a band above 3 kHz, which means that the speech signal can be isolated with a standard low-pass filter (low-pass filter).

Development of a digital filter

We will develop a digital filter with suitable characteristics. To create a low-pass filter, we use the function digitalfilter in which we specify the type of filter response lowpass with a transmission frequency of 3000 Hz, and we also announce a method for its synthesis. Elliptic this corresponds to an elliptical filter with an infinite impulse response (IIR) and the required characteristics (filter order - 10, ripple in the passband - no more than 0.1 dB, suppression in the boom band - no less than 60 dB):

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

Calculate and display the filter characteristics. First of all, we are interested in the amplitude-frequency response (AFC), which allows us to understand which frequencies of the spectrum of the original signal our filter will pass and which will suppress. To do this, use the function 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]:

Now we will display the phase-frequency response (FCH) of the filter with the function phaseresp:

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

Finally, we will calculate and display the impulse response of (THEIR) filter (more precisely, its first 100 samples, since our filter has an infinite number of THEM) by the function impresp:

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

Let's apply the developed filter to the original audio signal. filt and we will display the result of processing in the frequency domain.:

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]:

On the spectrum, you can verify that frequencies above 3 kHz have been successfully suppressed.

We will also display the result of filtering in the time domain.:

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

It can be seen from the graph that there is no noise filling between the spoken words. Let's verify the success of the processing by listening to the filtered audio signal.:

In [ ]: 
audioplayer(filtered, fs)

The applied filter successfully removed noise from the signal, but also suppressed the high-frequency components of the voice signal, which is why it does not sound so brightly colored with overtones and hissing when listening.

Conclusion

Using the example of the selected audio file, we examined the typical workflow of digital signal processing in Engee, which includes import, visualization, analysis, and processing.