Цифровая обработка сигналов в Engee на примере типового проекта с исходными данными

В данном примере мы разберём типовой проект по обработке, фильтрации и анализу данных в Engee. Данные будут загружены из бинарных файлов, определяющих подстилающую поверхность бетона и характеристики выходного сигнала.

Для начала обработки выполним чтение значения амплитуд исходного сигнала из .ana-файла:

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

Pkg.add("DelimitedFiles")

   Resolving package versions...
  No Changes to `~/.project/Project.toml`
  No Changes to `~/.project/Manifest.toml`

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]

Теперь считаем значения частоты входного сигнала из .xml-файла:

Pkg.add("XMLDict")

   Resolving package versions...
  No Changes to `~/.project/Project.toml`
  No Changes to `~/.project/Manifest.toml`

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

Опираясь на полученное значение частоты сигнала, создаем вектор временных отсчетов.

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

0.0:0.0004:31.6996

Выполним сохранение рассчитанных значений и параметров сигнала в MAT-файл.

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)

Загружаем коэффициенты заранее созданного фильтра низких частот:

# Загрузка данных из файла 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

Рассмотрим условия, при которых нам необходимо изменить частоту дискретизации входного сигнала с 2500 Гц до 200 Гц. Найдем соотношение этих частот, чтобы задать кратность изменения частоты дискретизации и сформировать корректный алгоритм фильтрации.

rat = rationalize(200/2500) 

up = rat.num
down = rat.den

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

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

Частота дискретизации сперва увеличивается в 2 раза, а затем уменьшается в 25 раз. Выполним фильтрацию при понижении частоты дискретизации.

После фильтрации рассчитаем новые временные отсчёты и построим результирующий график сравнения амплитуд входного и выходного сигнала.

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)

Как мы видим, длина во времени выходного сигнала значительно меньше, при этом общие тенденции поведения сигнала сохранены.

Теперь построим классический Фурье-спектр для проведения частотного анализа входного и выходного сигнала.

# Расчёт спектра сигнала
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)

Как мы можем видеть анализируя спектр, выходной сигнал имеет меньшую амплитуду и более узкий диапазон частот, чем входной сигнал.

Далее выполним вейвлет-преобразование. Оно используется для анализа временных рядов, таких как сигналы, звуки или изображения, с целью выявления их локальных характеристик во времени и частоте.

Оно позволяет представить сигнал одновременно в двух доменах – временном и частотном, что делает его особенно полезным при анализе нестационарных сигналов, где частота изменяется со временем.

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

using ContinuousWavelets, Wavelets

   Resolving package versions...
  No Changes to `~/.project/Project.toml`
  No Changes to `~/.project/Manifest.toml`
   Resolving package versions...
  No Changes to `~/.project/Project.toml`
  No Changes to `~/.project/Manifest.toml`

Разберём детально код, описанный ниже.

wavelet — это функция, которая создает объект волнового преобразователя. В нашем случае она принимает три аргумента.

1. Morlet — волновая функция Морле с параметром π. Функция Морле широко применяется в анализе сигналов благодаря своей способности эффективно выделять особенности сигнала на разных масштабах.
2. averagingType=NoAve() определяет тип усреднения. Здесь выбрано отсутствие усреднения. Это означает, что каждая точка в результирующем массиве будет соответствовать значению волнового коэффициента без дополнительного сглаживания.
3. β — параметр, определяющий ширину окна функции Морле. Чем выше значение β, тем более узкое окно задано и тем точнее локализация сигнала во временной области, но хуже разрешение в частотной области.

cwt — эта функция выполняет непрерывное волновое преобразование (CWT) исходного сигнала с использованием созданного волнового преобразователя.

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

Выполним частотно-временное преобразование Фурье. В Engee под такие преобразования нет стандартных функций, и их необходимо реализовывать вручную. Вот подробная схема, как это делается.

# Окно Хэмминга
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 = 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))

Опираясь на графики выше, мы можем определить:

1. какие частоты преобладают в нашем сигнале;
2. как эти частоты распределены относительно времени.

Теперь выполним анализ кепстра.

Кепстр (cepstrum) — это метод анализа сигналов, основанный на преобразовании Фурье логарифмического спектра сигнала. Кепстр широко применяется в обработке речевых сигналов, акустике и анализе вибрационных процессов. Существует несколько видов кепстров, каждый из которых предоставляет свою уникальную информацию о сигнале.

В Engee для выполнения анализа кепстра можно использовать пакеты FFTW и стандартные функции для выполнения преобразования Фурье. Далее представлены реализованные нами функции.

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)

Информация, предоставляемая комплексным кепстром:

1. Комплексный кепстр помогает разделить источники сигнала, особенно когда они имеют различные временные задержки. Это полезно, например, в системах эхо-подавления или идентификации многолучевого распространения.
2. Поскольку комплексный кепстр включает в себя информацию о фазе, он может использоваться для выявления и коррекции фазовых искажений в сигналах.

Информация, предоставляемая обратным комплексным кепстром, следующая.

1. Благодаря сохранению информации о фазе обратный комплексный кепстр позволяет точно восстанавливать исходный сигнал после обработки.
2. Обратный комплексный кепстр полезен для удаления шумов и помех, поскольку он сохраняет при этом важные компоненты сигнала.

Информация, предоставляемая действительным кепстром, такова.

1. Действительный кепстр часто используется для оценки основного тона речи или музыки, поскольку пики в кепстре соответствуют периодическим компонентам сигнала.
2. Действительный кепстр может эффективно подавлять аддитивный белый шум, так как он сглаживает логарифмический спектр.

Вывод

В данном примере мы раскрыли принципы, используемые для цифровой обработки сигналов, и применили различные методы анализы этих сигналов для решения разных инженерных задач.