Digital signal processing in Engee on the example of a typical project with initial data

In this example, we will walk through a typical data processing, filtering and analysis project in Engee. The data will be loaded from binary files defining the underlying concrete surface and the characteristics of the output signal.

To start processing, we will read the amplitude values of the original signal from the .ana file:

Pkg.add(["XMLDict", "ContinuousWavelets", "DelimitedFiles", "DSP", "MAT", "Wavelets"])

Pkg.add("DelimitedFiles")

using DelimitedFiles

filename = joinpath(("$(@__DIR__)/signal/sig0001.ana"))
io = open(filename, "r") # Открываем файл в режиме чтения двоичных данных
A = reinterpret(Float32, read(io)) # Читаем весь файл целиком в массив типа Float32
close(io) # Закрываем файл

println("Количество элементов: ", length(A))
println("Первые несколько элементов массива:")
show(A[1:min(length(A), 5)])

Количество элементов: 79250
Первые несколько элементов массива:
Float32[0.17833632, 0.091151185, 0.0041207145, 0.13451618, -0.12640864]

Now read the frequency values of the input signal from the .xml file:

Pkg.add("XMLDict")

using XMLDict

# Открываем файл
filename = joinpath("$(@__DIR__)/signal/sig0001.xml")  # Предполагаю, что nom - строковая переменная с путём к директории
xml_data = read(filename, String)

# Находим индекс начала строки с частотой
k = findfirst("frequency", xml_data)
if isnothing(k)
    error("Строка 'frequency' не найдена.")
end

# Извлекаем значение частоты
i = k.start + 11
ss = Char[]
while xml_data[i] != '"'
    push!(ss, xml_data[i])
    i += 1
end

# Преобразуем строку в число
fd = parse(Float64, join(ss))

# Выводим результат
println("Частота: ", fd)

Частота: 2500.0

Based on the obtained signal frequency value, we create a vector of time samples.

t = 0:(1/fd):((1/fd)*(length(A)-1))

0.0:0.0004:31.6996

Let's save the calculated values and signal parameters to MAT-file.

using MAT

# Преобразуйте t в Vector
vectorized_t = collect(t)

# Создайте словарь с данными
variables = Dict(
    "A" => A,
    "fd" => fd,
    "vectorized_t" => vectorized_t  
)

# Сохранить данные в файл .mat
filename = joinpath(@__DIR__, "signal/res.mat")
matwrite(filename, variables)

Load the coefficients of the previously created low-pass filter:

# Загрузка данных из файла f2500.mat
file = matopen(joinpath(@__DIR__, "signal", "f2500.mat"))  # Используем joinpath для безопасного построения пути
vars = read(file)
close(file)

# Доступ ко всем загруженным данным
# for var in vars
#     println(var)
# end

# Получение значения Num2
Num2 = vec(vars["Num2"])

3098-element Vector{Float64}:
 -0.00135077207137034
  0.00030632310930407626
  0.00027098304458891245
  0.00023859837725492806
  0.00020806857273815513
  0.00018030781760094202
  0.00015420602134399715
  0.00013069560414379414
  0.00010865108520715402
  8.903221313386013e-5
  7.068592854163531e-5
  5.461377664609944e-5
  3.961954370077718e-5
  ⋮
  5.461377664609944e-5
  7.068592854163531e-5
  8.903221313386013e-5
  0.00010865108520715402
  0.00013069560414379414
  0.00015420602134399715
  0.00018030781760094202
  0.00020806857273815513
  0.00023859837725492806
  0.00027098304458891245
  0.00030632310930407626
 -0.00135077207137034

Consider the conditions under which we need to change the sampling frequency of the input signal from 2500 Hz to 200 Hz. Let's find the ratio of these frequencies to set the multiplicity of sampling frequency change and form a correct filtering algorithm.

rat = rationalize(200/2500) 

up = rat.num
down = rat.den

print("Числитель и знаменатель: $rat")

Числитель и знаменатель: 2//25

The sampling frequency is first increased by a factor of 2 and then decreased by a factor of 25. Let's perform filtering at decreasing sampling frequency.

After filtering, calculate the new time samples and plot the resulting graph comparing the amplitudes of the input and output signal.

using DSP
using FFTW

new_a = vcat([A zeros(length(A))]...) # Создаем новый массив new_a, чередуя элементы A с нулями
diskr = DSP.conv(new_a, Num2) # Фильтрация

new_A = [diskr[1]] # Уменьшаем частоту дискретизации
for i in 1+down:down:length(diskr)
    push!(new_A, diskr[i])
end

new_t = 0:(down/(up*fd)):(length(new_A)-1)*(down/(up*fd))

plot(t,A)
plot!(new_t,new_A)

As we can see, the length in time of the output signal is much smaller, while the general tendencies of the signal behaviour are preserved.

Now let us construct a classical Fourier spectrum to perform frequency analysis of the input and output signal.

# Расчёт спектра сигнала
function compute_spectrum(signal, fs)
    n = length(signal)
    spectrum = abs.(fft(signal)) / n
    freqs = (0:n-1) .* (fs / n)
    spectrum[1:Int(n/2)], freqs[1:Int(n/2)]  # Вернуть половину спектра (для удобства)
end

compute_spectrum (generic function with 1 method)

spectrum_A, freqs_A = compute_spectrum(A, 2500)
plot(freqs_A, spectrum_A, xlabel="Frequency (Hz)", ylabel="Amplitude", legend=:topright)

spectrum_A, freqs_A = compute_spectrum(new_A, 200)
plot(freqs_A, spectrum_A, xlabel="Frequency (Hz)", ylabel="Amplitude", legend=:topright)

As we can see by analysing the spectrum, the output signal has smaller amplitude and narrower frequency range than the input signal.

Next, let's perform the wavelet transform. It is used to analyse time series such as signals, sounds or images in order to identify their local characteristics in time and frequency.

It allows the signal to be represented simultaneously in two domains, time and frequency, which makes it particularly useful in analysing non-stationary signals where frequency varies with time.

Pkg.add("Wavelets")
Pkg.add("ContinuousWavelets")

using ContinuousWavelets, Wavelets

Let's analyse the code described below in detail.

wavelet is a function that creates a wave converter object. In our case it takes three arguments.

1. morlet is a Morlet wave function with the parameter π. The Morlet function is widely used in signal analysis due to its ability to effectively extract signal features at different scales.
2. averagingType=NoAve() defines the type of averaging. No averaging is selected here. This means that each point in the resulting array will correspond to the value of the wave coefficient without additional smoothing.
3. β is a parameter determining the width of the Morlet function window. The higher the value of β, the narrower the window is set and the more accurate the signal localisation in the time domain, but the worse the resolution in the frequency domain.

cwt - this function performs continuous wave transformation (CWT) of the source signal using the created wave converter.

c = wavelet(Morlet(π*2), averagingType=NoAve(), β=3)
res = cwt(new_A, c)
heatmap(abs.(res)', xlabel= "time",ylabel="frequency")
xlims!(0, 3000)

Let's perform the Fourier frequency-time transform. There are no standard functions in Engee for such transformations, and they must be implemented manually. Here is a detailed scheme of how to do it.

# Окно Хэмминга
function hamming(N::Int)
    n = 0:N-1
    return 0.54 .- 0.46 .* cos.(2π * n / (N - 1))
end

hamming (generic function with 1 method)

l=6000 # Размер окна
window=vcat(hamming(l),zeros(length(new_A)-l))

P = periodogram(new_A; fs=200, window=window, nfft=length(new_A))
Power_log = 20*log10.(power(P));
Freq = freq(P);

# Визуализация спектра
plot(Freq, Power_log, xlabel="Frequency (Hz)", ylabel="Power Spectrum (dB)", title="Power Spectrum", linewidth=2)

# Вычисление спектограммы
S = spectrogram(new_A, fs=200,)

# Построение графика
heatmap(S.time, S.freq, S.power, xlabel="Время, с", ylabel="Частота, Гц", title="Диаграмма мощности сигнала", colorbar_title="Мощность, м²/c⁴",xlim =(0, 15))

Based on the graphs above, we can determine:

1. what frequencies dominate our signal;
2. how these frequencies are distributed with respect to time.

Now let's perform a cepstra analysis.

Cepstrum is a signal analysis method based on the Fourier transform of the logarithmic spectrum of a signal. Cepstrum is widely used in speech signal processing, acoustics and vibration analysis. There are several types of cepstrums, each providing its own unique information about the signal.

In Engee, FFTW packages and standard functions for performing the Fourier transform can be used to perform cepstral analysis. The following are the functions we have implemented.

using FFTW  # Для выполнения быстрого преобразования Фурье

# Комплексный кепстр
function cceps(x::Vector{Float64})
    X = fft(x,1)                     # Фурье-преобразование
    p = atan.(imag(X), real(X))
    log_mag = complex.(log.(abs.(X)),p) # Логарифм модуля спектра (добавляем малое значение для стабильности)
    cepstrum = real(ifft(log_mag,1)) # Обратное преобразование Фурье
    return cepstrum
end
# Обратный комплексный кепстр
function icceps(c::Vector{Float64})
    exp_spectrum = exp.(fft(c,1))    # Экспонента спектра
    signal = real(ifft(exp_spectrum,1)) # Обратное преобразование Фурье
    return signal
end
# Действительный кепстр
function rceps(x::Vector{Float64})
    X = fft(x,1)
    log_mag = log.(abs.(X) .+ 1e-10) # Логарифм модуля (с добавлением стабилизатора)
    cepstrum = real(ifft(log_mag)) 
    return cepstrum
end

cepstrum = cceps(new_A)
restored_signal = abs.(icceps(cepstrum))
real_cepstrum = abs.(rceps(new_A))

# # Визуализация
p1 = plot(abs.(cepstrum), title="Cepstrum", ylabel="Magnitude")
p2 = plot(abs.(restored_signal), title="Restored Signal", ylabel="Amplitude")
p3 = plot(abs.(real_cepstrum), title="Real Cepstrum", ylabel="Amplitude")
plot(p1, p2, p3, layout=(3,1), legend=false)

The information provided by the complex cepstrum:

1. The complex cepstrum helps to separate signal sources, especially when they have different time delays. This is useful, for example, in echo cancellation or multipath propagation identification systems.
2. Since the complex cepstrum includes phase information, it can be used to detect and correct phase distortion in signals.

The information provided by the inverse complex cepstrum is as follows. 1- By storing the phase information, the inverse complex cepstrum allows the original signal to be accurately reconstructed after processing. 2. Inverse complex cepstrum is useful for removing noise and interference because it preserves important signal components.

The information provided by a valid cepstrum is as follows.

1. The valid cepstrum is often used to estimate the fundamental tone of speech or music because the peaks in the cepstrum correspond to the periodic components of the signal.
2. The valid cepstrum can effectively suppress additive white noise because it smooths the logarithmic spectrum.

Conclusion

In this example, we have uncovered the principles used for digital signal processing and applied various methods of analysing these signals to solve different engineering problems.