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; # Library of graph rendering
using Noise; # Noise Addition Library
using DigitalComm; # Library of communication systems
using DSP; # Digital Signal Processing Library
using FFTW; # The Fourier Transform Library
plotlyjs(); # Enabling the rendering method
Setting the sine wave
All transformations in the example will be performed on this sinusoid.
fs = 5000; # Steps per second
t = [0:1/fs:1;]; # Time
y = sin.(2*pi*10*t); # Setting the sine wave
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="frequency, Hz", ylabel="spectral power density", 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); # Quantization of the original sinusoid
# Output of the result
plot(y[1:250], label="the original signal")
plot!(y_quant[1:250], label="quantized signal at 80 levels")
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="frequency, Hz", ylabel="spectral power density", 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="frequency, Hz", ylabel="spectral power density", 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 (that's why it's 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); # Adding noise
plot(t_new, signal, label="a sinusoid with noise")
plot!(t_new, y_d, label="a sine wave without noise")
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="the original signal")
plot!(Out_sig, label="filtered signal")
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="the original signal")
plot!(Out_sig2, label="filtered signal")
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 the assignment of 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) # Bandwidth detection
designmethod = FIRWindow(hanning(128)) # Determining the window size
Out_fir = filt(digitalfilter(responsetype, designmethod), Out_sig)
plot(Out_sig, label="the original signal")
plot!(Out_fir, label="the original signal")
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="frequency, Hz", ylabel="spectral power density", label="filt")
plot!(freq(y_out)[1:25], power(y_out)[1:25], xlabel="frequency, Hz", ylabel="spectral power density", label="conv")
plot!(freq(y_fir)[1:25], power(y_fir)[1:25], xlabel="frequency, Hz", ylabel="spectral power density", 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="frequency, Hz", ylabel="spectral power density", 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="frequency, Hz", ylabel="spectral power density", legend=false)
As we can see, the frequency of the filtered signal has halved.
println("The size of the entrance: ",size(Out_fir,1))
println("Downsampling size: ",size(rs,1))
println("Upsampling size: ",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.