应用 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 滤波器具有有限脉冲响应，这使得它们可以稳定求解，易于设计。它们没有内部反馈。

我们应用的滤波器类型是截止频率为 5 Hz 的低通滤波器（Lowpass）。它使用汉宁窗方法设计的长度为 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")

滤波后，信号变得平滑，噪声成分被抑制。现在的波形能更好地反映低频成分。

噪声发射峰值几乎消失，这证明高频成分的滤波效果良好。

32 赫兹以上的频率被完全抑制。频谱的主要振幅集中在 0 赫兹（直流分量）和 3 赫兹。

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

低通滤波器（低通滤波器）

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

我们使用的是截止频率为 0.2 Hz 的低通滤波器。它采用四阶椭圆滤波器方法设计，通带公差为 0.5 dB，阻带公差为 30 dB。
使用多项式表示法（PolynomialRatio）实现传递函数。

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

由于主信号中的残余噪声清晰可见，因此该滤波器的性能比前一个滤波器更差。虽然有平滑处理，但其质量并不令人满意，因为噪音仍然对最终结果产生重大影响。

处理后信号的频谱显示，在 130 Hz 以下的频率存在峰值，这表明对高频成分的过滤不足：滤波器未能有效抑制截止频率以上的频率。

椭圆滤波器具有较高的瞬态响应陡度，但在通带中会出现明显的纹波，本例中就出现了这种情况。

截止频率（0.2 Hz）和所选滤波器参数对该信号来说可能不是最佳选择。此外，滤波器阶数（4 阶）可能不足以达到所需的噪声抑制效果。

高通滤波器（高通滤波器）

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

我们在这里使用的滤波器类型是截止频率为 40 Hz 的高通滤波器，使用 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")

经过滤波后，信号与噪声原始信号几乎没有区别，这表明定义信号主要波动的低频信息已被完全去除。

频谱显示，40 赫兹以下的所有频率（包括 24 赫兹以下的重要低频成分）都被完全抑制。这导致主要信号信息丢失，只剩下高频噪声。

截止频率为 40 Hz 的高通滤波器完全不适合分析该信号，因为信号的重要频率（0.1 Hz）位于低频范围，而滤波器去除了低频。

带通滤波器（Band-pass filters）

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

采样频率为1000 Hz。

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

滤波器平滑了噪声，这在将处理后的信号与原始噪声信号进行比较时可以明显看出。然而，信号低频波动的有用信息却丢失了。只剩下信号高频部分的 "片段"，而这部分并不包含重要信息。

处理后信号的频谱显示，滤波器有效地抑制了 10 赫兹以下和 70 赫兹以上的频率，但保留了 10 至 40 赫兹的范围，在这种情况下，该范围并不包含有用的信息。这导致信号在这一频段被噪声成分所代表，失去了原有的低频特性。

带阻滤波器（Band-stop filters）

这些滤波器会衰减指定频率范围内的信号，允许所有其他频率的信号通过。

在本例中，我们使用一个四阶巴特沃斯滤波器来实现带阻滤波器，抑制 10-40 Hz 范围内的频率。

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 赫兹至 32 赫兹范围内的频率被抑制。包括噪声在内的其他频率成分几乎没有变化。

这就解释了为什么信号的噪声特征得以保留。

结论

我们研究了各种滤波器来解决这个问题，发现汉宁窗 FIR 滤波器是最佳选择。它能够平滑有效地去除特定截止频率以上的噪声，使我们能够分离出有用的低频信息（0.1 赫兹），而这正是我们的主要目标。

综合上述比较，我们可以得出以下结论：对于此类任务，值得使用截止频率与有用信息范围最匹配的低通滤波器。