以原始数据典型项目为例进行工程中的数字信号处理
在这个例子中，我们将分析一个典型的项目，用于处理，过滤和分析Engee中的数据。 数据将从定义底层混凝土表面和输出信号特征的二进制文件下载。
要开始处理，我们将从原始信号中读取振幅值。ana档案。:
In [ ]:
# Pkg.add(["XMLDict", "ContinuousWavelets", "DelimitedFiles", "DSP", "MAT", "Wavelets"])
In [ ]:
Pkg.add("DelimitedFiles")
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)])
现在我们正在计算从输入信号的频率值。xml文件:
In [ ]:
Pkg.add("XMLDict")
In [ ]:
using XMLDict
# 打开文件
filename = joinpath("$(@__DIR__)/signal/sig0001.xml") # 我假设nom是一个带有目录路径的字符串变量。
xml_data = read(filename, String)
# 我们用频率找到行开头的索引
k = findfirst("frequency", xml_data)
if isnothing(k)
error("没有找到字符串"频率"。")
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)
基于接收到的信号频率值，我们创建时间样本的向量。
In [ ]:
t = 0:(1/fd):((1/fd)*(length(A)-1))
Out[0]:
将信号的计算值和参数保存到MAT文件中。
In [ ]:
using MAT
# 将t转换为矢量
vectorized_t = collect(t)
# 使用数据创建字典
variables = Dict(
"A" => A,
"fd" => fd,
"vectorized_t" => vectorized_t
)
# 保存数据到a。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]:
让我们考虑需要将输入信号的采样频率从2500Hz更改为200Hz的条件。 我们将找到这些频率的比值，以便设置采样频率变化的频率并形成正确的滤波算法。
In [ ]:
rat = rationalize(200/2500)
up = rat.num
down = rat.den
print("分子和分母：$rat")
采样率首先增加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]:
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
我们来详细分析下面描述的代码。
wavelet是创建波形转换器对象的函数。 在我们的例子中，它需要三个参数。
- Morlet是参数π的Morlet波函数。 Morlet函数由于能够有效隔离不同尺度下的信号特征而被广泛应用于信号分析。
- averagingType=NoAve()定义求平均值的类型。 这里不选择平均值。 这意味着所得数组中的每个点将对应于波浪系数的值，而无需额外的平滑。
- β是确定Morlaix函数窗口宽度的参数。 Β值越高，窗口设置得越窄，信号在时域的定位越准确，但在频域的分辨率越差。
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]:
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="容量,m2/c",xlim =(0, 15))
Out[0]:
根据上面的图表，我们可以确定:
- 什么频率在我们的信号中占主导地位;
- 这些频率如何相对于时间分布。
现在让我们进行kepstra分析。
倒谱是一种基于信号对数谱傅立叶变换的信号分析方法。 Cepstr广泛应用于语音信号处理、声学和振动分析。 有几种类型的上限，每种类型都提供有关信号的独特信息。
在Engee中，用于执行傅立叶变换的FFTW包和标准函数可用于执行kepstra分析。 以下是我们实现的功能。
In [ ]:
using FFTW # 执行快速傅立叶变换
# 综合capstr
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
# 反向复数kepstr
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]:
集成kepstrom提供的信息:
- 集成的capstr有助于分离信号源，特别是当它们具有不同的时间延迟时。 这例如在回声消除或多径识别系统中是有用的。
- 由于复数capstr包含相位信息，因此可用于检测和校正信号中的相位失真。
逆复数kepstrom提供的信息如下。
- 由于相位信息的保存，反向复数倒谱允许您在处理后准确地恢复原始信号。
- 反向复数滤波器可用于消除噪声和干扰，因为它保留了重要的信号分量。
有效的kepstrom提供的信息如下。
- 有效的kepstr通常用于估计语音或音乐的音高，因为kepstr中的峰值对应于信号的周期性分量。
- 有效的倒谱可以有效抑制加性白噪声，因为它可以平滑对数谱。
结论
在这个例子中，我们揭示了用于数字信号处理的原理，并应用了分析这些信号的各种方法来解决各种工程问题。