Digital signal processing in Engee using the example of a typical project with raw data
In this example, we will analyze a typical project for processing, filtering and analyzing data in Engee. The data will be downloaded from binary files that define 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)])
Now we are counting 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)
Based on the received signal frequency value, we create a vector of time samples.
t = 0:(1/fd):((1/fd)*(length(A)-1))
Save the calculated values and parameters of the signal to a 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)
Loading the coefficients of a pre-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"])
Let's 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 in order to set the frequency of sampling frequency variation and form a correct filtering algorithm.
rat = rationalize(200/2500)
up = rat.num
down = rat.den
print("Числитель и знаменатель: $rat")
The sampling rate first increases by 2 times, and then decreases by 25 times. Let's perform filtering when the sampling rate is lowered.
After filtering, we will calculate new time samples and plot the resulting graph comparing the amplitudes of the input and output signals.
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 time length of the output signal is much shorter, while the general trends of the signal behavior are preserved.
Now we will construct a classical Fourier spectrum for frequency analysis of the input and output signals.
# Расчёт спектра сигнала
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
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 analyzing the spectrum, the output signal has a smaller amplitude and a narrower frequency range than the input signal.
Next, we perform the wavelet transform. It is used to analyze time series such as signals, sounds, or images in order to identify their local characteristics in time and frequency.
It allows you to represent a signal simultaneously in two domains – time and frequency, which makes it especially useful in the analysis of non-stationary signals, where the frequency varies with time.
Pkg.add("Wavelets")
Pkg.add("ContinuousWavelets")
using ContinuousWavelets, Wavelets
Let's analyze in detail the code described below.
wavelet is a function that creates a waveform converter object. In our case, it takes three arguments.
- Morlet is the Morlet wave function with the parameter π. The Morlet function is widely used in signal analysis due to its ability to effectively isolate signal features at different scales.
- 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.
- β is a parameter that determines the width of the Morlaix function window. The higher the β value, the narrower the window is set and the more accurate the localization of the signal in the time domain, but the worse the resolution in the frequency domain.
CWT — This function performs continuous wave transformation (CWT) of the original signal using a built-in 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 time-frequency Fourier transform. There are no standard functions in Engee for such transformations, and they must be implemented manually. Here is a detailed diagram of how it is done.
# Окно Хэмминга
function hamming(N::Int)
n = 0:N-1
return 0.54 .- 0.46 .* cos.(2π * n / (N - 1))
end
l=6000 # Размер окна
window=vcat(hamming(l),zeros(length(new_A)-l))
P = DSP.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 = DSP.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:
- which frequencies prevail in our signal;
- How these frequencies are distributed relative to time.
Now let's perform the kepstra analysis.
Cepstrum is a signal analysis method based on the Fourier transform of the logarithmic spectrum of a signal. Cepstr is widely used in speech signal processing, acoustics, and vibration analysis. There are several types of caps, each of which provides its own unique information about the signal.
In Engee, FFTW packages and standard functions for performing the Fourier transform can be used to perform kepstra analysis. The following are the functions implemented by us.
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)
Information provided by the integrated kepstrom:
- An integrated capstr helps to separate the signal sources, especially when they have different time delays. This is useful, for example, in echo cancellation or multipath identification systems.
- Since the complex capstr includes phase information, it can be used to detect and correct phase distortions in signals.
The information provided by the inverse complex kepstrom is as follows.
- Due to the preservation of phase information, the reverse complex cepstr allows you to accurately restore the original signal after processing.
- The reverse complex filter is useful for removing noise and interference, as it preserves important signal components.
The information provided by a valid kepstrom is as follows.
- A valid kepstr is often used to estimate the pitch of speech or music, since the peaks in the kepstr correspond to periodic components of the signal.
- A valid cepstr can effectively suppress additive white noise, as it smooths the logarithmic spectrum.
Conclusion
In this example, we have revealed the principles used for digital signal processing and applied various methods of analyzing these signals to solve various engineering problems.