应用DSP库过滤器

为了演示来自DSP库的不同类型的过滤器，我们准备了每个过滤器的简要描述和应用示例。

在这里，我们来看看基本的过滤器。:

具有有限脉冲响应的滤波器,
低通滤波器,
高频滤波器。

您可以在此示例中详细了解构建过滤器的所有原则：https://engee.com/community/ru/catalogs/projects/sravnitelnyi-analiz-filtrov

让我们从连接这个库开始。值得注意的是，DSP库区分滤波器系数和有状态滤波器。

Pkg.add("DSP")
using DSP
using FFTW

┌ Warning: could not download http://pkgserver.pkgserver.svc.cluster.local/registries
│   exception = RequestError: HTTP/1.1 404 Not Found while requesting http://pkgserver.pkgserver.svc.cluster.local/registries
└ @ Pkg.Registry /usr/local/julia/share/julia/stdlib/v1.10/Pkg/src/Registry/Registry.jl:69
   Resolving package versions...
  No Changes to `~/.project/Project.toml`
  No Changes to `~/.project/Manifest.toml`

我们将宣布输入信号，稍后我们将用它来测试我们的过滤器。

我们还将声明一个辅助函数。

Fs = 1000  # Частота дискретизации
t = 0:1/Fs:1-1/Fs  # Временной интервал
x = sin.(2π .* 0.1 .* t) .+ randn(length(t))  # Сигнал с шумом (0.1 Гц компонент и шум)
plot(x)

# Расчёт спектра сигнала
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_x, freqs_x = compute_spectrum(x, Fs)
plot(freqs_x, spectrum_x, xlabel="Frequency (Hz)", ylabel="Amplitude", legend=:topright)

接下来，让我们来看看各种过滤器的实现示例。

FIR滤波器（有限脉冲响应）

FIR滤波器具有有限的脉冲响应，这使得它们可以稳定求解并且易于设计。他们没有内部反馈。

我们使用的滤波器类型是截止频率为5Hz的低通滤波器。它使用由汉宁窗法设计的长度为64的FIR滤波器应用于信号。

responsetype = Lowpass(5)
designmethod = FIRWindow(hanning(64))
y = filt(digitalfilter(responsetype, designmethod; fs=Fs), x)

plot(x, label="Original Signal", title="Signal")
plot!(y, label="Filtered Signal", xlabel="Sample", ylabel="Amplitude")

spectrum_y, freqs_y = compute_spectrum(y, Fs)

plot(freqs_x, spectrum_x, label="Original Signal Spectrum", xlabel="Frequency (Hz)", ylabel="Amplitude", legend=:topright)
plot!(freqs_y, spectrum_y, label="Filtered Signal Spectrum", title="Signal Spectrum")

经过滤波后，信号变得平滑，抑制了噪声成分。波形现在更好地反映了低频分量。

峰值噪声辐射实际上已经消失，这证实了高频分量的良好滤波。

高于32Hz的频率被完全抑制。频谱的主要振幅集中在0Hz（直流分量）和3Hz的频率。

这表明滤波器有效地抑制了截止频率以上的噪声。

低通滤波器

这些滤波器通过频率低于设定阈值的信号，并衰减频率较高的信号。

我们使用的滤波器类型是截止频率为0.2hz的低通滤波器。采用四阶椭圆滤波器方法设计，通带公差偏差为0.5db，势垒带公差偏差为30db。
它被实现为使用多项式表示（多项式比率）的传递函数。

responsetype = Lowpass(0.2)
designmethod = Elliptic(4, 0.5, 30)
tf = convert(PolynomialRatio, digitalfilter(responsetype, designmethod))
numerator_coefs = coefb(tf)
denominator_coefs = coefa(tf)

y = filt(tf, x)

plot(x, label="Original Signal", title="Signal")
plot!(y, label="Filtered Signal", xlabel="Sample", ylabel="Amplitude")

spectrum_y, freqs_y = compute_spectrum(y, Fs)

plot(freqs_x, spectrum_x, label="Original Signal Spectrum", xlabel="Frequency (Hz)", ylabel="Amplitude", legend=:topright)
plot!(freqs_y, spectrum_y, label="Filtered Signal Spectrum", title="Signal Spectrum")

由于残留噪声在主信号中可见，因此该滤波器比前一滤波器更能应对这一任务。有平滑，但其质量不令人满意，因为干扰继续对最终结果产生显着影响。

处理信号的频谱显示在高达130Hz的频率处存在峰值，这表明高频分量的滤波不足：滤波器无效地抑制高于其截止频率的频率。

椭圆滤波器具有高陡度的过渡特性，但是它们可以遭受带宽中的显着波纹，这在这种情况下被观察到。

截止频率(0.2Hz)和所选滤波器参数可能不是该信号的最佳选择。此外，滤波器阶数(第4)可能不足以实现所需的干扰抑制。

高通滤波器

这些滤波器通过频率高于设定阈值的信号，并衰减频率较低的信号。

我们在此应用的滤波器类型是截止频率为40Hz的高通滤波器，使用4阶巴特沃斯滤波器实现。

responsetype = Highpass(40)  # Высокочастотный фильтр, отсечка на 40 Гц, fs = 1000 Гц
designmethod = Butterworth(4)         # Баттерворт, порядок 4
highpass_filter = digitalfilter(responsetype, designmethod; fs=Fs)

y = filt(highpass_filter, x)  # Применяем фильтр к сигналу x

plot(x, label="Original Signal", title="Signal")
plot!(y, label="Filtered Signal", xlabel="Sample", ylabel="Amplitude")

spectrum_y, freqs_y = compute_spectrum(y, Fs)

plot(freqs_x, spectrum_x, label="Original Signal Spectrum", xlabel="Frequency (Hz)", ylabel="Amplitude", legend=:topright)
plot!(freqs_y, spectrum_y, label="Filtered Signal Spectrum", title="Signal Spectrum")

经过滤波后，信号实际上与有噪声的原始信号没有区别，这表明确定信号主要波动的低频信息已被完全删除。

频谱显示，所有高达40Hz的频率（包括高达24Hz的显着低频分量）都被完全抑制。这导致有关信号的基本信息丢失，只留下高频噪声。

截止频率为40Hz的高通滤波器完全不适合分析该信号，因为信号的显着频率(0.1Hz)位于滤波器消除的低频范围内。

带通滤波器

这些滤波器通过给定频带内的信号，衰减该频带以外的信号。
我们在这里使用的滤波器类型是带宽为10至40Hz的带通滤波器，使用4阶巴特沃斯滤波器实现。

采样率：1000Hz。

responsetype = Bandpass(10, 40)
designmethod = Butterworth(4)
y = filt(digitalfilter(responsetype, designmethod; fs=Fs), x)

plot(x, label="Original Signal", title="Signal")
plot!(y, label="Filtered Signal", xlabel="Sample", ylabel="Amplitude")

spectrum_y, freqs_y = compute_spectrum(y, Fs)

plot(freqs_x, spectrum_x, label="Original Signal Spectrum", xlabel="Frequency (Hz)", ylabel="Amplitude", legend=:topright)
plot!(freqs_y, spectrum_y, label="Filtered Signal Spectrum", title="Signal Spectrum")

滤波器已经平滑了噪声，这在将处理后的信号与原始噪声信号进行比较时是显而易见的。然而，关于信号的低频波动的有用信息已经丢失。信号的高频部分只有一个"片段"，它不携带显着的信息。

处理信号的频谱表明，滤波器有效地抑制了低于10Hz和高于70Hz的频率，但保留了从10到40Hz的范围，在这种情况下不包含有用的信息。这导致信号由该频带中的噪声分量表示，失去了其原始的低频特性。

带阻滤波器

这些滤波器衰减给定频率范围内的信号，跳过所有其他频率。

在这种情况下，我们正在实施一个抑制10-40Hz范围内频率的带通滤波器，并使用4阶巴特沃斯滤波器实现。

responsetype = Bandstop(10, 40)  # Полосо-заградительный фильтр, подавляет частоты от 10 до 40 Гц, fs = 1000 Гц
designmethod = Butterworth(4)            # Баттерворт, порядок 4

y = filt(digitalfilter(responsetype, designmethod; fs=Fs), x)  # Применяем фильтр к сигналу x (определите x заранее)


plot(x, label="Original Signal", title="Signal")
plot!(y, label="Filtered Signal", xlabel="Sample", ylabel="Amplitude")

spectrum_y, freqs_y = compute_spectrum(y, Fs)

plot(freqs_x, spectrum_x, label="Original Signal Spectrum", xlabel="Frequency (Hz)", ylabel="Amplitude", legend=:topright)
plot!(freqs_y, spectrum_y, label="Filtered Signal Spectrum", title="Signal Spectrum")

滤波后，信号看起来与原始噪声信号几乎相同，这表明滤波器对信号可见结构的轻微影响。

频谱显示在12至32Hz范围内的频率抑制。其余的频率成分，包括噪声，几乎保持不变。

这解释了信号的持续噪声性质。

结论

我们研究了各种滤波器来解决这个问题，发现带有汉宁窗口的FIR滤波器是最好的选择。它能够平稳有效地消除给定截止频率以上的噪声，从而隔离有用的低频信息（0.1Hz），这是主要任务。

总结您所做的比较，我们可以得出以下结论：对于此类任务，值得使用截止频率与有用信息范围最佳对应的低通滤波器。