自适应50Hz滤波器

Question

我应用以下滤波器来去除信号上的50Hz净噪声：

#python code
def filter_50(signal):
    for i in np.arange(50,500,50):
        fs = 1000.0  # Sample frequency (Hz)
        f0 = i  # Frequency to be removed from signal (Hz)
        w0 = f0 / (fs / 2)  # Normalized Frequency
        Q= 30
        b, a = iirnotch(w0, Q)
        signal = scipy.signal.filtfilt(b, a, signal)
    return(signal)

我用Q玩过，但它不够好。

Answer 1

嗯，使用滤波器，你总是在信号失真和去除不需要的频率之间折衷。滤波后总会有某种信号剩余，这取决于滤波器衰减系数。如果指定为陷波滤波器，巴特沃思滤波器可以有几乎100%的衰减。下面是使用巴特沃斯滤波器的效果：

这显示的原始信号是50赫兹，目标是如果滤波器足够好，我们应该看不到任何信号滤波后。然而，在使用带宽为15 Hz的二阶巴特沃斯滤波器后，我们确实看到了一些信号，尤其是在信号的开始和结束时，这是由于滤波器失真造成的。

并且滤波器的频率响应在频域（幅度和相位）看起来是这样的。

因此，虽然相位变化平滑，但巴特沃斯滤波器幅度的“陷波”效应也是平滑的。

改变Q将改变50赫兹的衰减水平和失真。我认为总的来说，如果您的噪声接近或重叠感兴趣的信号，使用iirnotch是一个好主意，否则Butterwoth可能是一个更好的选择。

以下是这些数字的代码：

from scipy.signal import filtfilt, iirnotch, freqz, butter
from scipy.fftpack import fft, fftshift, fftfreq
import numpy as np
from matplotlib import pyplot


def do_fft(y, fs):
    Y = fftshift(fft(y, 2 ** 12))
    f = fftshift(fftfreq(2 ** 12, 1 / fs))
    return f, Y


def make_signal(fs, f0, T=250e-3):
    # T is total signal time
    t = np.arange(0, T, 1 / fs)
    y = np.sin(2 * np.pi * f0 * t)
    return t, y


def make_plot():
    fig, ax = pyplot.subplots(1, 2)
    ax[0].plot(t, y)
    ax[0].plot(t, y_filt)
    ax[0].set_title('Time domain')
    ax[0].set_xlabel('time [seconds]')
    ax[1].plot(f, abs(Y))
    ax[1].plot(f, abs(Y_filt))
    ax[1].set_title('Frequency domain')
    ax[1].set_xlabel('Freq [Hz]')

    # filter response
    fig, ax = pyplot.subplots(1, 2)
    ax[0].plot(filt_freq, abs(h))
    ax[0].set_title('Amplitude')
    ax[0].set_xlim([0, 200])
    ax[0].set_xlabel('Freq [Hz]')
    ax[1].plot(filt_freq, np.unwrap(np.angle(h)) * 180 / np.pi)
    ax[1].set_title('Phase')
    ax[1].set_xlim([0, 200])
    ax[1].set_xlabel('Freq [Hz]')

    pyplot.show()


fs = 1000
f0 = 50
t, y = make_signal(fs=fs, f0=f0)
f, Y = do_fft(y, fs=1000)

# Filtering using iirnotch
w0 = f0/(fs/2)
Q = 30
b, a = iirnotch(w0, Q)

# filter response
w, h = freqz(b, a)
filt_freq = w*fs/(2*np.pi)
y_filt = filtfilt(b, a, y)
f, Y_filt = do_fft(y_filt, fs)

make_plot()


w0 = [(f0-15)/(fs/2), (f0+15)/(fs/2)]
b, a = butter(2, w0, btype='bandstop')
w, h = freqz(b, a)
filt_freq = w*fs/(2*np.pi)
y_filt = filtfilt(b, a, y)
f, Y_filt = do_fft(y_filt, fs)
make_plot()