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.
We will need basic libraries for digital signal processing and working with WAV files.:
neededLibs = ["DSP", "WAV", "Base64"]
for lib in neededLibs
try
eval(Meta.parse("using $lib"))
catch ex
Pkg.add(lib)
eval(Meta.parse("using $lib"))
end
end
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.
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.:
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
Let's listen to the original audio file (a useful signal against a background of high-frequency noise):
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).
plot(t[1:10:end], original[1:10:end],
xlabel = "Time, from",
ylabel = "The amplitude",
title = "The original signal in time")
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
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 = "The spectrum of the original signal")
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
spgr = DSP.spectrogram(original, 511, 256; fs = fs);
heatmap(spgr.time, spgr.freq, pow2db.(spgr.power),
ylim=(0,10000),
clim=(-140,-40),
xlabel = "Time, seconds",
ylabel = "Frequency, Hz",
title = "Spectrogram of the original signal")
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 (LPF).
We will develop a digital filter with suitable characteristics. To create a low-pass filter, we use the function digitalfilterlowpassElliptic
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
H, w = freqresp(myfilt);
freq_vec = fs*w/(2*pi);
plot(freq_vec, pow2db.(abs.(H))*2, linewidth=3,
xlabel = "Frequency, Hz",
ylabel = "Amplitude, dB",
title = "Frequency response of the low-pass filter")
Now we will display the phase-frequency response (FCH) of the filter with the function phaseresp
phi, w = phaseresp(myfilt);
plot(freq_vec, rad2deg.(phi/2), linewidth = 3,
xlabel = "Frequency, Hz",
ylabel = "Phase, degree",
title = "Low-pass filter frequency response")
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
myimpresp = impresp(myfilt,100);
plot(myimpresp, line=:stem, marker=:circle,
xlabel = "THEIR counts",
ylabel = "THEIR values",
title = "Low-frequency impulse response")
Let's apply the developed filter to the original audio signal. 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 = "Frequency, Hz",
ylabel = "Мощность/частота, дБ/Гц",
title = "The spectrum of the processed signal")
On the spectrum, it can be verified that frequencies above 3 kHz have been successfully suppressed.
We will also display the filtering result in the time domain.:
plot(t[1:10:end], [original[1:10:end], filtered[1:10:end]],
xlabel = "Time, from",
ylabel = "The amplitude",
title = "Processed signal in time")
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.:
audioplayer(filtered, fs)
The applied filter successfully removed noise from the signal, but also suppressed the high-frequency components of the voice signal, which makes it sound less brightly colored with overtones and hissing when listening.
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.