Engee documentation
Notebook

Filtering in Engee

In this example, we will look at the possibilities of specifying signals, methods of filtering them and methods of changing the signal frequency, as well as analyse the estimation of the power spectral density of signals at each of the stages of interaction with the signal.

Connecting libraries

In [ ]: 
Pkg.add(["Noise", "DSP", "DigitalComm"])

   Resolving package versions...
    Updating `~/.project/Project.toml`
  [a7b11256] + DigitalComm v1.2.0
  [81d43f40] + Noise v0.3.3
    Updating `~/.project/Manifest.toml`
  [a7b11256] + DigitalComm v1.2.0
  [81d43f40] + Noise v0.3.3
  [d96e819e] + Parameters v0.12.3
  [e409e4f3] + PoissonRandom v0.4.4
  [3a884ed6] + UnPack v1.0.2
Precompiling project...
  ✓ Noise
  1 dependency successfully precompiled in 4 seconds. 23 already precompiled.
In [ ]: 
using Plots; # Билиотека отрисовки графиков
using Noise; # Библиотека добавления шумов
using DigitalComm; # Библиотека систем связи
using DSP; # Библиотека цифровой обработки сигналов
using FFTW; # Библиотека преобразований Фурье
plotlyjs(); # Подключения метода отрисовки

Setting a sine wave

All the transformations in the example will be performed on this sinusoid.

In [ ]: 
fs = 5000; # Шагов в секунду
t = [0:1/fs:1;]; # Время

y = sin.(2*pi*10*t); # Задание синусоиды

plot(t,y)

┌ Warning: attempting to remove probably stale pidfile
│   path = "/user/.jlassetregistry.lock"
└ @ Pidfile /usr/local/ijulia-core/packages/Pidfile/DDu3M/src/Pidfile.jl:260
Out[0]:

Let's estimate the power spectral density of the signal through periodogram - this is a function of periodogram construction, based on the calculation of the square of the modulus of the Fourier transform for the data sequence.

Power spectral density (PSD) in physics and signal processing is a function describing the distribution of signal power as a function of frequency, i.e. the power attributable to a unit frequency interval.

In [ ]: 
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)
Out[0]:

Quantisation

The next stage of our signal conversion is quantisation - dividing the range of reference values of the signal into a finite number of levels and rounding these values to one of the two levels nearest to them.

In [ ]: 
y_quant = quantization(y, 80, minv=-1, maxv=1); # Квантование исходной синусоиды

# Вывод результата
plot(y[1:250], label="исходный сигнал")
plot!(y_quant[1:250], label="квантованный сигнал по 80 уровням")
Out[0]:
In [ ]: 
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)
Out[0]:

As we can see, this transformation does not affect the power spectral density of the signal and does not shift its frequency.

Data thinning

The next action we will perform on the original sine wave is data thinning. We will take every second sample from the original signal.

In [ ]: 
y_d = y_quant[1:2:end]
t_new = t[1:2:end]
plot(t_new,y_d)
Out[0]:
In [ ]: 
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)
Out[0]:

The periodogram shows that data thinning increases the frequency of the signal. This is due to the fact that due to the reduction of the number of samples by a factor of two, the frequency of oscillations has correspondingly increased by a factor of two.

Noise superposition

Noise is unwanted phenomena that prevent us from getting information from a useful signal.

Noise is random in nature.

In digital signal processing, additive white Gaussian noise (ABGN) is most commonly used.

Its characteristics are:

  • uniform power spectral density (that is why it is called white);
  • normal distribution of time values (hence it is called Gaussian);
  • summed with the useful signal (therefore it is called additive);
  • ABGS is statistically independent of the useful signal.
In [ ]: 
signal = addNoise(y_d, 20); # Добавление шумов

plot(t_new, signal, label="синусоида с шумом")
plot!(t_new, y_d, label="синусоида без шума")
Out[0]:

It is important to note that the function addNoise defines a tuple from which we need only the first vector. The following code shows the data types of input addNoise, its output and the result of extracting the first vector from the tuple.

In [ ]: 
typeof(y_d)
Out[0]:
Vector{Float64} (alias for Array{Float64, 1})
In [ ]: 
typeof(signal)
Out[0]:
Tuple{Vector{Float64}, Vector{Float64}}
In [ ]: 
signal=signal[1];
In [ ]: 
typeof(signal)
Out[0]:
Vector{Float64} (alias for Array{Float64, 1})

Signal filtering in Engee

In this example we will look at three filters, but you can find any filter you need in the Engee documentation. The first method we will consider is to specify the filter through the filter coefficients, in the form of two vectors a and b, then realise the filter representation through the coefficients of the numerator b and denominator a of the 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 highest degree to lowest degree.

In [ ]: 
b = [0.75, 0.5];
a = [1, 0.2, 0.1];
f = PolynomialRatio(b,a);

Out_sig = filt(f,signal);

This filter is not a decimating filter. This can be verified by comparing the size of the filter input and output.

In [ ]: 
size(Out_sig,1)==size(signal,1)
Out[0]:
true
In [ ]: 
plot(signal, label="исходный сигнал")
plot!(Out_sig, label="отфильтрованный сигнал")
Out[0]:

The second variant of specifying the same filter is through the impulse response using the function impresp, and then calculating the filter output using the function conv.

In [ ]: 
h = impresp(f);
Out_sig2 = conv(signal,h);

plot(signal, label="исходный сигнал")
plot!(Out_sig2, label="отфильтрованный сигнал")
Out[0]:

Let's compare the output of the function filt and the output of the function conv.

In [ ]: 
plot(Out_sig[1:2500]-Out_sig2[1:2500])
Out[0]:

As we can see from the graph, the outputs of the function filt and the output of the function conv almost completely coincide.

The third filtering variant that we will consider is setting a window FIR filter. It will be different from the previous ones, because it has a delay relative to the input signal equal to the window size.

In [ ]: 
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="исходный сигнал")
Out[0]:

To compare the filter outputs we will use the Welch periodogram.

In [ ]: 
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")
Out[0]:

As we see, the main differences between the filter outputs lie in their power spectral density.

Changing the frequency of the signal

In order to increase or decrease the number of samples, Engee uses the command resample, which allows you to execute upsampling and downsampling of the signal.

Note: upsampling and downsampling are specified through the relation n//k, where:

n: the number of points to be created;

k: the number of nearest neighbours to be considered for each point.

In [ ]: 
# 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)
Out[0]:

As we can see, the frequency of the filtered signal has doubled.

In [ ]: 
# 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)
Out[0]:

As we can see, the frequency of the filtered signal has halved.

In [ ]: 
println("Размер входа: ",size(Out_fir,1))
println("Размер при downsampling: ",size(rs,1))
println("Размер при upsampling: ",size(rs2,1))

Размер входа: 2501
Размер при downsampling: 1251
Размер при upsampling: 5001
In [ ]: 
plot(rs)
plot!(rs2)
Out[0]:

The same effects can be observed when comparing signals and their sizes. At downsampling the number of samples of the signal decreases, and at upsampling, on the contrary, it increases relative to the original signal.

Conclusion

In this example we have demonstrated the possibilities of interaction with signals and showed that Engee can be fully applicable to digital signal processing.