在 0.667 和 3.833 Hz 范围内快速实现 IIR 滤波器答案

【问题标题】：Swift implementation of IIR filter in the 0.667 and 3.833 Hz range在 0.667 和 3.833 Hz 范围内快速实现 IIR 滤波器
【发布时间】：2018-05-06 21:42:29
【问题描述】：

我需要去除 0.667 - 3.833 Hz 范围之外的频率（心率范围：40 到 230 bpm）。我查看了 Apple 文档，但我只看到了 AVFoundation Butterworth 过滤器。显然，“值的有效范围是 20 Hz 到 (SampleRate/2)。”该滤波器来自 AVAudioUnitEQ，因此可以理解...我似乎无法找到在我正在使用的频率范围内工作的任何滤波器。

我是否必须实现自己的 Butterworth 过滤器，或者我缺少什么？

编辑：所以我对文档进行了深入研究，发现了 IIR filter 和 FIR filter，但这并不是我所期望的...... IIR 滤波器是“级联双二阶 IIR 滤波器”，不是Butterworth filter（这是我更喜欢的）。我做了一个little reading，但我不知道这个过滤器应该如何工作。如果我使用这些系数，我可以实现与巴特沃斯滤波器类似的功能吗？如果是这样，我该怎么做？目标是获得没有波纹的平坦通带。

我正在关注的博文：http://www.ignaciomellado.es/blog/Measuring-heart-rate-with-a-smartphone-camera

包含一些双二阶滤波器建议的博文：http://www.ignaciomellado.es/blog/Measuring-heart-rate-with-a-smartphone-camera

【问题讨论】：

如果您使用 Swift for iOS 或 macOS 进行编码，Accelerate 框架 (vDSP) 包含双二阶 IIR 过滤例程。但是，一个数值更稳定的解决方案将涉及首先进行低通滤波，然后在带通滤波之前对数据进行下采样。
@hotpaw2 是的，这正是我正在做的。我刚刚阅读了一篇提到您的建议的博客文章。很有可能，我会这样做。
@JacquesThibodeau 您是否能够在 Swift 中找到可接受的解决方案？

标签： swift signal-processing

【解决方案1】：

IIR 过滤器的实现在任何语言中都非常简单明了。 IIR 的设计是挑战和有趣的部分。可以进行很多权衡。双二阶 IIR 滤波器是非常好的二阶滤波器（可以使用巴特沃斯系数进行配置），对于许多应用来说，一个就足够了，但是您也可以级联它们以制作更高阶的滤波器，如果您需要。

首先要弄清楚你真正需要什么，这通常意味着你需要了解你想要过滤的噪声在频域中是什么样的，然后设计一个抑制噪声的滤波器。

周围有很多工具可以让您了解给定 IIR 系数集的频率响应是什么样的。

这里是one，来自快速谷歌搜索。

一个关键参数是您对数据源进行采样的频率（采样率）

一旦有了系数，实现就非常简单，它只是当前输入样本和许多先前输入和输出的函数。

你可以找到一些细节here，但这里有一个函数。

y[n] = b0* x[n] + b1 * x[n-1] + b2 * x[n-2] - a1 * y[n-1] - a2 * y[n-2]

其中 x[n] 是您当前的样本，x[n-1] 是前一个样本等，y[n] 是您当前的输出值，而 y[n-1], y[n- 2] 是以前的输出。

如果您真的想深入了解过滤器设计的细节，我发现的最佳资源之一是使用SciPy，同样有很多关于如何做到这一点的教程。

我没有看过它，但这里有一个 youtube video 描述如何在 Python 中设计巴特沃斯过滤器。

正如我之前所说，一旦有了系数，任何语言的实现都非常简单。例如，这是我在 c++ 中做的一个非常简单的 Biquad，用于过滤 Arduino 项目中的数据。这可以很容易地转换为 swift 或任何其他语言。

class BiQuad {

public:
  BiQuad(double b0, double b1, double b2, double a1, double a2)
      : b0_(b0), b1_(b1), b2_(b2), a1_(a1), a2_(a2) {}


  double processSample(double xn) {
    xn_ = xn;
    yn_ = b0_ * xn_ + b1_ * xnm1_ + b2_ * xnm2_ - a1_ * ynm1_ - a2_ * ynm2_;
    ynm2_ = ynm1_;
    ynm1_ = yn_;
    xnm2_ = xnm1_;
    xnm1_ = xn_;
    return yn_;
  }

  double getCurrentY() { return yn_; }

  void prime(double value) {
    xn_ = xnm1_ = xnm2_ = yn_ = ynm1_ = ynm2_ = value;
  }

private:
  double b0_;
  double b1_;
  double b2_;
  double a1_;
  double a2_;
  double xn_;
  double xnm1_ = 0.0;
  double xnm2_ = 0.0;
  double yn_ = 0.0;
  double ynm1_ = 0.0;
  double ynm2_ = 0.0;
};

希望这会有所帮助。

【讨论】：