Engee 文档
Notebook

以带有初始数据的典型项目为例，介绍 Engee 中的数字信号处理技术

在本示例中，我们将通过 Engee 演示一个典型的数据处理、滤波和分析项目。数据将从二进制文件中加载，这些文件定义了底层混凝土表面和输出信号的特征。

开始处理时，我们将从 .ana 文件中读取原始信号的振幅值：

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

   Resolving package versions...
  No Changes to `~/.project/Project.toml`
  No Changes to `~/.project/Manifest.toml`
In [ ]: 
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 文件中读取输入信号的频率值：

In [ ]: 
Pkg.add("XMLDict")

   Resolving package versions...
  No Changes to `~/.project/Project.toml`
  No Changes to `~/.project/Manifest.toml`
In [ ]: 
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

根据获得的信号频率值，我们创建一个时间采样向量。

In [ ]: 
t = 0:(1/fd):((1/fd)*(length(A)-1))
Out[0]:
0.0:0.0004:31.6996

将计算值和信号参数保存到 MAT 文件中。

In [ ]: 
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)

加载之前创建的低通滤波器的系数：

In [ ]: 
# Загрузка данных из файла 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"])
Out[0]:
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 Hz 变为 200 Hz。让我们找出这些频率的比值，以设定采样频率变化的倍数，并形成正确的滤波算法。

In [ ]: 
rat = rationalize(200/2500) 

up = rat.num
down = rat.den

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

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

首先将采样频率提高 2 倍，然后降低 25 倍。让我们以递减的采样频率进行滤波。

滤波后，计算新的时间采样，并绘制输入和输出信号振幅比较图。

In [ ]: 
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)
Out[0]:

我们可以看到，输出信号的时间长度大大缩短，而信号行为的总体趋势却得以保留。

现在，让我们构建一个经典的傅里叶频谱，对输入和输出信号进行频率分析。

In [ ]: 
# Расчёт спектра сигнала
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
Out[0]:
compute_spectrum (generic function with 1 method)
In [ ]: 
spectrum_A, freqs_A = compute_spectrum(A, 2500)
plot(freqs_A, spectrum_A, xlabel="Frequency (Hz)", ylabel="Amplitude", legend=:topright)
Out[0]:
In [ ]: 
spectrum_A, freqs_A = compute_spectrum(new_A, 200)
plot(freqs_A, spectrum_A, xlabel="Frequency (Hz)", ylabel="Amplitude", legend=:topright)
Out[0]:

通过频谱分析我们可以看到，输出信号的振幅比输入信号小，频率范围比输入信号窄。

接下来，让我们进行小波变换。它用于分析信号、声音或图像等时间序列，以识别它们在时间和频率上的局部特征。

它允许同时在时间和频率两个域中表示信号，因此在分析频率随时间变化的非稳态信号时特别有用。

In [ ]: 
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`

让我们来详细分析一下下面描述的代码。

小波是一个创建波形转换器对象的函数。在我们的例子中，它需要三个参数。

  1. morlet 是参数为 π 的莫列特波函数。Morlet 函数能有效提取不同尺度的信号特征，因此被广泛应用于信号分析。 averagingType=NoAve() 定义了平均化类型。此处选择无平均。这意味着结果数组中的每个点都将对应于波形系数的值，而无需额外的平滑处理。 β 是决定莫雷特函数窗口宽度的参数。β 值越大，窗口越窄，时域信号定位越精确，但频域分辨率越差。

cwt - 该函数使用创建的波形转换器对源信号执行连续波变换（CWT）。

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

让我们执行傅立叶频率-时间变换。Engee 中没有用于此类变换的标准函数，必须手动执行。下面是详细的实现方法。

In [ ]: 
# Окно Хэмминга
function hamming(N::Int)
    n = 0:N-1
    return 0.54 .- 0.46 .* cos.(2π * n / (N - 1))
end
Out[0]:
hamming (generic function with 1 method)
In [ ]: 
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)
Out[0]:
In [ ]: 
# Вычисление спектограммы
S = DSP.spectrogram(new_A, fs=200,)

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

根据上述图表，我们可以确定

  1. 哪些频率主导我们的信号；
  2. 这些频率是如何随时间分布的。

现在我们来进行 cepstra 分析。

倒频谱是一种基于信号对数频谱傅里叶变换的信号分析方法。倒频谱广泛应用于语音信号处理、声学和振动分析。倒频谱有多种类型，每种类型都能提供有关信号的独特信息。

在 Engee 中，FFTW 软件包和执行傅立叶变换的标准函数可用于进行倒频谱分析。以下是我们已经实现的函数。

In [ ]: 
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)
Out[0]:

复倒频谱提供的信息： 1.复倒谱有助于分离信号源，尤其是当信号源具有不同的时间延迟时。例如，这在回声消除或多径传播识别系统中非常有用。 2.由于复倒谱包含相位信息，因此可用于检测和纠正信号中的相位失真。

反复倒频谱提供的信息如下。 1- 通过存储相位信息，逆复倒频谱可以在处理后准确地重建原始信号。 2 反复倒频谱由于保留了重要的信号成分，因此有助于去除噪声和干扰。

有效倒频谱提供的信息如下。 1.有效倒频谱通常用于估计语音或音乐的基音，因为倒频谱中的峰值与信号的周期成分相对应。 2.有效倒频谱能有效抑制加性白噪声，因为它能平滑对数频谱。

结论

在这个例子中，我们揭示了数字信号处理的原理，并应用了各种分析这些信号的方法来解决不同的工程问题。