Filtering in Engee
In this example, we will analyze the possibilities of specifying signals, ways to filter them, and methods for changing the frequency of the signal, as well as analyze the estimation of the spectral power density of the signals at each stage of interaction with the signal.
Connecting libraries
Pkg.add(["Noise", "DSP", "DigitalComm"])
using Plots; # Билиотека отрисовки графиков
using Noise; # Библиотека добавления шумов
using DigitalComm; # Библиотека систем связи
using DSP; # Библиотека цифровой обработки сигналов
using FFTW; # Библиотека преобразований Фурье
plotlyjs(); # Подключения метода отрисовки
Setting the sine wave
All transformations in the example will be performed on this sinusoid.
fs = 5000; # Шагов в секунду
t = [0:1/fs:1;]; # Время
y = sin.(2*pi*10*t); # Задание синусоиды
plot(t,y)
Let's estimate the spectral power density of the signal using periodogram, which is a periodogram construction function based on calculating the square of the Fourier transform modulus for a sequence of data.
The spectral power density (SPM) in physics and signal processing is a function describing the distribution of signal power depending on frequency, that is, the power per unit frequency interval.
y_p = DSP.periodogram(y, onesided=true, nfft=length(y), fs=fs);
plot(freq(y_p)[1:20], power(y_p)[1:20], xlabel="частота, Гц", ylabel="спектральная плотность мощности", legend=false)
Quantization
The next step in converting our signal will be its quantization — splitting the range of reference values of the signal into a finite number of levels and rounding these values to one of the two levels closest to them.
y_quant = quantization(y, 80, minv=-1, maxv=1); # Квантование исходной синусоиды
# Вывод результата
plot(y[1:250], label="исходный сигнал")
plot!(y_quant[1:250], label="квантованный сигнал по 80 уровням")
y_p = DSP.periodogram(y_quant, onesided=true, nfft=length(y_quant), fs=fs);
plot(freq(y_p)[1:20], power(y_p)[1:20], xlabel="частота, Гц", ylabel="спектральная плотность мощности", legend=false)
As we can see, this transformation does not affect the spectral power density of the signal and does not shift its frequency.
Thinning of data
The next step that we will perform with the original sinusoid is data thinning. We will take every second sample from the original signal.
y_d = y_quant[1:2:end]
t_new = t[1:2:end]
plot(t_new,y_d)
y_p = DSP.periodogram(y_d, onesided=true, nfft=length(y_d), fs=fs);
plot(freq(y_p)[1:20], power(y_p)[1:20], xlabel="частота, Гц", ylabel="спектральная плотность мощности", legend=false)
The periodogram shows that thinning the data increases the frequency of the signal. This is due to the fact that due to the halving of the number of samples, the oscillation frequency has doubled accordingly.
Noise overlay
Noise is an undesirable phenomenon that prevents us from receiving information from a useful signal.
Noise is random in nature.
In digital signal processing, additive white Gaussian noise (ABGSH) is most often used.
Its characteristics are:
- Uniform spectral power density (which is why it is called white);
- the normal distribution of time values (therefore it is called Gaussian);
- it is combined with a useful signal (therefore it is called additive);
- The ABSH is statistically independent of the useful signal.
signal = addNoise(y_d, 20); # Добавление шумов
plot(t_new, signal, label="синусоида с шумом")
plot!(t_new, y_d, label="синусоида без шума")
It is important to note that the function addNoise sets a tuple from which we need only the first vector.
The following code shows the types of login data addNoise, its output and the result of extracting the first vector from the tuple.
typeof(y_d)
typeof(signal)
signal=signal[1];
typeof(signal)
Signal filtering in Engee
In this example, we will analyze three filters, but in the Engee documentation you can find any filter you need.
The first method that we will consider is setting the filter through the filter coefficients, in the form of two vectors. a and b, then we implement the representation of the filter in terms of the coefficients of the numerator
b and the denominator a transfer function using the function PolynomialRatio contained in the library DSP.
Note:
b and a can be specified as objects Polynomial, or
vectors ordered from the largest degree to the smallest.
b = [0.75, 0.5];
a = [1, 0.2, 0.1];
f = PolynomialRatio(b,a);
Out_sig = filt(f,signal);
This filter is not decemming. You can verify this by comparing the sizes of the filter input and output.
size(Out_sig,1)==size(signal,1)
plot(signal, label="исходный сигнал")
plot!(Out_sig, label="отфильтрованный сигнал")
The second option for setting the same filter is through the impulse response using the function impresp, and then calculating the filter output using the function conv.
h = impresp(f);
Out_sig2 = conv(signal,h);
plot(signal, label="исходный сигнал")
plot!(Out_sig2, label="отфильтрованный сигнал")
Let's compare the output of the filt function and the output of the conv function.
plot(Out_sig[1:2500]-Out_sig2[1:2500])
As we can see from the graph, the outputs of the function are filt and the function output conv they almost completely coincide.
The third filtering option that we will analyze is setting a window FIR filter. It will differ from the previous ones, as it has a delay relative to the input signal equal to the window size.
responsetype = Lowpass(5; fs=30) # Определение полосы пропускания
designmethod = FIRWindow(hanning(128)) # Определение размеров окна
Out_fir = filt(digitalfilter(responsetype, designmethod), Out_sig)
plot(Out_sig, label="исходный сигнал")
plot!(Out_fir, label="исходный сигнал")
To compare the filter outputs, we use the Welch periodogram.
n=div(length(Out_sig),8);
y_out = welch_pgram(Out_sig, n, div(n,2), onesided=true, nfft=length(Out_sig), fs=fs);
n=div(length(Out_sig2),8);
y_in = welch_pgram(Out_sig2, n, div(n,2), onesided=true, nfft=length(Out_sig2), fs=fs);
n=div(length(Out_fir),8);
y_fir = welch_pgram(Out_fir, n, div(n,2), onesided=true, nfft=length(Out_fir), fs=fs);
plot(freq(y_in)[1:25], power(y_in)[1:25], xlabel="частота, Гц", ylabel="спектральная плотность мощности", label="filt")
plot!(freq(y_out)[1:25], power(y_out)[1:25], xlabel="частота, Гц", ylabel="спектральная плотность мощности", label="conv")
plot!(freq(y_fir)[1:25], power(y_fir)[1:25], xlabel="частота, Гц", ylabel="спектральная плотность мощности", label="fir")
As we can see, the main differences between the filter outputs lie in their spectral power density.
Changing the signal frequency
In order to increase or decrease the number of counts, Engee uses the command resample, which allows you to perform upsampling and downsampling the signal.
Note: upsampling and downsampling they are set using the relation n//k, where:
n: the number of points that need to be created;
k: the number of nearest neighbors to consider for each point.
# Downsampling
rs = resample(Out_fir, 1//2)
n=div(length(rs),8);
y_out = welch_pgram(rs, n, div(n,2), onesided=true, nfft=length(rs), fs=fs);
plot(freq(y_out)[1:25], power(y_out)[1:25], xlabel="частота, Гц", ylabel="спектральная плотность мощности", legend=false)
As we can see, the frequency of the filtered signal has doubled.
# Upsampling
rs2 = resample(Out_fir, 2//1)
n=div(length(rs2),8);
y_out = welch_pgram(rs2, n, div(n,2), onesided=true, nfft=length(rs2), fs=fs);
plot(freq(y_out)[1:25], power(y_out)[1:25], xlabel="частота, Гц", ylabel="спектральная плотность мощности", legend=false)
As we can see, the frequency of the filtered signal has halved.
println("Размер входа: ",size(Out_fir,1))
println("Размер при downsampling: ",size(rs,1))
println("Размер при upsampling: ",size(rs2,1))
plot(rs)
plot!(rs2)
We can observe the same effects when comparing the signals and their sizes. By downsampling the number of signal samples decreases, and when upsampling on the contrary, it increases relative to the original signal.
Conclusion
In this example, we demonstrated the possibilities of interacting with signals and showed that Engee can be fully applied to digital signal processing.