以原始数据典型项目为例进行工程中的数字信号处理

在这个例子中，我们将分析一个典型的项目，用于处理，过滤和分析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

让我们考虑需要将输入信号的采样频率从2500Hz更改为200Hz的条件。让我们找到这些频率的比率，以便设置采样频率变化的频率并形成正确的滤波算法。

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是创建波形转换器对象的函数。在我们的例子中，它需要三个参数。

Morlet是参数π的Morlet波函数。 Morlet函数由于能够有效隔离不同尺度下的信号特征而被广泛应用于信号分析。
averagingType=NoAve()定义求平均值的类型。这里不选择平均值。这意味着所得数组中的每个点将对应于波浪系数的值，而无需额外的平滑。
β是确定Morlaix函数窗口宽度的参数。 Β值越高，窗口设置得越窄，信号在时域的定位越准确，但在频域的分辨率越差。

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))

根据上面的图表，我们可以确定:

什么频率在我们的信号中占主导地位;
这些频率如何相对于时间分布。

现在让我们进行kepstra分析。

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

在Engee中，用于执行傅立叶变换的FFTW包和标准函数可用于执行kepstra分析。以下是我们实现的功能。

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)

集成kepstrom提供的信息:

集成的capstr有助于分离信号源，特别是当它们具有不同的时间延迟时。这例如在回声消除或多径识别系统中是有用的。
由于复数capstr包含相位信息，因此可用于检测和校正信号中的相位失真。

逆复数kepstrom提供的信息如下。

由于相位信息的保存，反向复数倒谱允许您在处理后准确地恢复原始信号。
反向复数滤波器可用于消除噪声和干扰，因为它保留了重要的信号分量。

有效的kepstrom提供的信息如下。

有效的kepstr通常用于估计语音或音乐的音高，因为kepstr中的峰值对应于信号的周期性分量。
有效的倒谱可以有效抑制加性白噪声，因为它可以平滑对数谱。

结论

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