Цифровая обработка сигналов в Engee на примере типового проекта с исходными данными
В данном примере мы разберём типовой проект по обработке, фильтрации и анализу данных в Engee. Данные будут загружены из бинарных файлов, определяющих подстилающую поверхность бетона и характеристики выходного сигнала.
Для начала обработки выполним чтение значения амплитуд исходного сигнала из .ana-файла:
# 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)])
Теперь считаем значения частоты входного сигнала из .xml-файла:
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)
Опираясь на полученное значение частоты сигнала, создаем вектор временных отсчетов.
t = 0:(1/fd):((1/fd)*(length(A)-1))
Выполним сохранение рассчитанных значений и параметров сигнала в 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"])
Рассмотрим условия, при которых нам необходимо изменить частоту дискретизации входного сигнала с 2500 Гц до 200 Гц. Найдем соотношение этих частот, чтобы задать кратность изменения частоты дискретизации и сформировать корректный алгоритм фильтрации.
rat = rationalize(200/2500)
up = rat.num
down = rat.den
print("Числитель и знаменатель: $rat")
Частота дискретизации сперва увеличивается в 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
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
Разберём детально код, описанный ниже.
wavelet — это функция, которая создает объект волнового преобразователя. В нашем случае она принимает три аргумента.
- Morlet — волновая функция Морле с параметром π. Функция Морле широко применяется в анализе сигналов благодаря своей способности эффективно выделять особенности сигнала на разных масштабах.
- averagingType=NoAve() определяет тип усреднения. Здесь выбрано отсутствие усреднения. Это означает, что каждая точка в результирующем массиве будет соответствовать значению волнового коэффициента без дополнительного сглаживания.
- β — параметр, определяющий ширину окна функции Морле. Чем выше значение β, тем более узкое окно задано и тем точнее локализация сигнала во временной области, но хуже разрешение в частотной области.
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
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))
Опираясь на графики выше, мы можем определить:
- какие частоты преобладают в нашем сигнале;
- как эти частоты распределены относительно времени.
Теперь выполним анализ кепстра.
Кепстр (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)
Информация, предоставляемая комплексным кепстром:
- Комплексный кепстр помогает разделить источники сигнала, особенно когда они имеют различные временные задержки. Это полезно, например, в системах эхо-подавления или идентификации многолучевого распространения.
- Поскольку комплексный кепстр включает в себя информацию о фазе, он может использоваться для выявления и коррекции фазовых искажений в сигналах.
Информация, предоставляемая обратным комплексным кепстром, следующая.
- Благодаря сохранению информации о фазе обратный комплексный кепстр позволяет точно восстанавливать исходный сигнал после обработки.
- Обратный комплексный кепстр полезен для удаления шумов и помех, поскольку он сохраняет при этом важные компоненты сигнала.
Информация, предоставляемая действительным кепстром, такова.
- Действительный кепстр часто используется для оценки основного тона речи или музыки, поскольку пики в кепстре соответствуют периодическим компонентам сигнала.
- Действительный кепстр может эффективно подавлять аддитивный белый шум, так как он сглаживает логарифмический спектр.
Вывод
В данном примере мы раскрыли принципы, используемые для цифровой обработки сигналов, и применили различные методы анализы этих сигналов для решения разных инженерных задач.