RC低通滤波器

转自：RC低通滤波器 - 知乎

RC low pass filter，或者叫直流滤波器(DC Filter)， 其作用是采用软件的方式，把一输入交流信号中的直流成分滤除掉。

其传递函数为:

其中K为加权因子，值域为 [0 ~ 1]，值越小由收敛越快，但信号的强度也会被大打折扣。通常该取接近于1的值，经常值为0.999。

下图为一组测试实例，图中黑色曲线为4组不同频率及振幅的正弦波的叠加，绿色曲线为经过直流滤波后的效果，其因子K=0.999；可以看到原始数据中的直波成分逐渐被过滤掉了。

上图制作的python代码：

import numpy as np
import matplotlib.pyplot as plt
import ctypes
import os, sys

os.chdir(os.path.abspath(os.path.dirname(sys.argv[0])))

'''
 Begin: C function and parameters and return declaration
'''
cLibs = ctypes.CDLL("{0}\\cLib\\cLib.dll".format(os.getcwd()))
cLibs.DoLowPass_Filter.restype = ctypes.c_float  # declare return type;
cLibs.DoHighPass_Filter.restype = ctypes.c_float  # declare return type;
cLibs.DoDc_Filter.restype = ctypes.c_float  # declare return type;
cLibs.DoSquared_Filter.restype = ctypes.c_float

'''
call init function.
'''
cLibs.InitDcFilter(ctypes.c_float(0.999))
cLibs.InitHighPass(ctypes.c_float(100),ctypes.c_float(16000))
cLibs.InitSquaredFilter(ctypes.c_float(0.99875))


def sin_wave(a, f, fs, phi, t, dc=0):
    """

    Generator sine wave sequence.
    若时间序列长度为 t=1s, 采样频率fs=1000Hz,则采样时间间隔 Ts=1/fs=0.001s;
    对于时间序列采样点个数为 n=t/Ts=1/0.001s=1000,即有1000个采样点，每个点的时间间隔为Ts.
    :param a: 振幅
    :param f: 信号频率
    :param fs: 采样频率
    :param phi: 相位
    :param t: 时间长度
    :param dc: 直流偏置
    :return: wave array.
    """
    ts = 1 / fs
    n = t / ts
    n = np.arange(n)
    y = dc + a * np.sin(2 * np.pi * f * n * ts + phi * (np.pi / 180))
    return y


'''
dc 滤波器
'''


def dc_filter(x_array, k=0.999):
    """

    :param x_array: 输入数据序列
    :param k: 滤波系数
    :return: 滤波后数据序列
    """
    y = []
    x_prev = 0
    y_prev = 0

    for x_val in x_array:
        y_prev = x_val - x_prev + k * y_prev
        x_prev = x_val
        y.append(y_prev)
    return y


'''
图形显示
'''
tspan = 0.5
Fs = 16000

hz_10 = sin_wave(a=0.3, f=10, fs=Fs, phi=0, t=tspan, dc=4)
hz_30 = sin_wave(a=0.3, f=30, fs=Fs, phi=0, t=tspan, dc=3)
hz_50 = sin_wave(a=0.3, f=50, fs=Fs, phi=0, t=tspan, dc=2)
hz_250 = sin_wave(a=0.5, f=250, fs=Fs, phi=0, t=tspan, dc=1)

hz_raw = np.array(hz_10) + np.array(hz_30) + np.array(hz_50) + np.array(hz_250)

hz_dc = []
hz_iir = []
hz_squared = []
hz_rooted = []
for itm in hz_raw:
    flted = cLibs.DoDc_Filter(ctypes.c_float(itm))
    hz_dc.append(flted)
    flted = cLibs.DoHighPass_Filter(ctypes.c_float(flted))
    hz_iir.append(flted)
    flted = cLibs.DoSquared_Filter(ctypes.c_float(flted))
    hz_squared.append(flted)
    flted = np.sqrt(flted)
    hz_rooted.append(flted)

x = np.arange(0, tspan, 1 / Fs)
plt.xlabel('t/s')
plt.ylabel('y')
plt.grid()
plt.plot(x, hz_raw, 'k')
plt.plot(x, hz_dc, linewidth=1, color='green')
plt.plot(x, hz_iir, 'r--')
plt.plot(x, hz_squared, 'b-.')
plt.plot(x, hz_rooted, 'red')
plt.legend(['raw', 'dc', 'iir', 'squared_exp_avg', 'squared_rooted'], loc=1, fontsize=26)
plt.show()