python-滤波

发表于 2020-11-26 更新于 2021-12-11 分类于 Language 阅读次数：本文字数： 6.7k 阅读时长 ≈ 6 分钟

香农采样定理是这样描述的：采样频率fs至少为被采集信号最高频率的2倍。少于 2 倍会出现混叠

采样频率少于2倍的信号频率时，会导致原本的高频信号被采样成低频信号，如下图所示。红色信号是原始的高频信号，但是由于采样频率不满足采样定理的要求，导致实际采样点如图中蓝色实点表示，将这些蓝色点连成曲线，可以明显的看出这是一个低频信号。在图示的时间长度内，红色信号有18个周期，但采样后的蓝色信号只有2个周期。也就是采样后的信号频率成分为原始信号频率成分的1/9。

这就是所谓的混叠。对连续信号进行等间隔采样时，如果采样频率不满足采样定理，采样后信号的频率就会发生混叠，即高于奈奎斯特频率的频率成分将被重构成低于奈奎斯特频率的信号。这种频谱的重叠导致的失真称为混叠，也就是高频信号被混叠成了低频信号。

滤波器存在滤波陡度，在滤波截止频率（奈奎斯特频率）以上的一些区域还存在混叠的可能性，这个区域对应的带宽的80%以上部分，也就是带宽的80%-100%区域。如下图所示，高于奈奎斯特频率以上的频率成分会关心奈奎斯特频率镜像到带宽的80%-100%区域，形成混叠，而带宽80%以内的区域，是无混叠的。

当按采样定理设置采样频率时，带宽的80%以上频带还可能存在混叠，如下图红框所示区域即遭受了频率混叠的影响。

既然采样定理要求的是2倍，那为什么要用2.56倍呢？基于以下两个方面的原因。

要关心的频带内无混叠

采样频率 fs ≥2.5fmax

如要求100Hz内无混叠，则采样频率应设置成250Hz，带宽为125Hz，带宽的80%为100Hz，因此，存在混叠可能性的带宽80%以上区域已位于感兴趣的频带之外了。当采样频率高于关心的最高频率2.5倍时，关心的频带内已无混叠了。

要方便计算机处理

快速傅立叶变换要求处理的数据块包含的数据点为2^N，而计算机也只能用0和1来存储数据，因此，计算机处理的数据时，如果是2^N会更方便些。我们知道256=2^8，因此，离2.5最近的2.56便成为了一个重要的“优先数”（先借用一下优先数这个概念）。

改变信号频率使输出信号降至最大值的0.707倍（对应-3dB），或0.5倍时（对应-6dB），该频率称为截止频率。 cut-off frequency
对于一个滤器器来说，在高频端和低频端各有一个截止频率，分别称为上截止频率fH和下截止频率fL。两个截止频率之间的频率范围称为通频带。 BW=fH-fL。

Chebyshev filter
切比雪夫滤波器

也被称为等涟波滤波器（equal ripple filter），是在通带或阻带上频率响应幅度等波纹波动的滤波器。在通带波动的为“I型切比雪夫滤波器”，在阻带波动的为“II型切比雪夫滤波器”。

Butterworth filter

scipy.signal.butter(N, Wn, btype=‘low’, analog=False, output=‘ba’)

输入参数：

N:滤波器的阶数
Wn：归一化截止频率。计算公式Wn=2*截止频率/采样频率。（注意：根据采样定理，采样频率要大于两倍的信号本身最大的频率，才能还原信号。截止频率一定小于信号本身最大的频率，所以Wn一定在0和1之间）。当构造带通滤波器或者带阻滤波器时，Wn为长度为2的列表。
btype : 滤波器类型{‘lowpass’, ‘highpass’, ‘bandpass’, ‘bandstop’},
output : 输出类型{‘ba’, ‘zpk’, ‘sos’},

输出参数：
b，a: IIR滤波器的分子（b）和分母（a）多项式系数向量。output=‘ba’
z,p,k: IIR滤波器传递函数的零点、极点和系统增益. output= ‘zpk’
sos: IIR滤波器的二阶截面表示。output= ‘sos’

N为滤波器阶数，wc为3dB截止频率，巴特沃斯低通滤波器的特点为：整个频带内单调下降，且非常光滑，并且阶数越高，阻带内下降的越陡峭、越干脆。无论阶数多高，都在经过（wc，3dB)一点。

Bessel filter

DE>0.707: Low amplitude distortion and high phase distortion (Chebyshev filter).
DE=0.707: Low amplitude distortion and moderate phase distortion (Butterworth filter).
DE=0.577: High altitude distortion and low phase distortion (Bessel filter)

巴特沃斯滤波器通带最平坦，阻带下降慢。
切比雪夫滤波器通带等纹波，阻带下降较快。
贝塞尔滤波器通带等纹波，阻带下降慢。也就是说幅频特性的选频特性最差。但是，贝塞尔滤波器具有最佳的线性相位特性。
此外，还有椭圆滤波器，椭圆滤波器在通带等纹波（阻带平坦或等纹波），阻带下降最快。

按照最佳逼近特性或者滤波通带特性分类，主要为巴特沃斯滤波器（Butterworth）、切比雪夫滤波器（Chebyshev）、贝塞尔滤波器（Bessel）和椭圆滤波器（Elliptic）四种

通过式滤波器可以让参考频率一侧的频率成分完全通过该滤波器，同时对另一侧的频率成分做线性的衰减，就是，一边让通过，一边逐渐被滤除。在信号学中，通过的区域被称为通带，滤除的区域被叫做阻带，在通过式滤波器中，参考频率通常被称为截止频率。

scipy.signal.butter(N, Wn, btype=’low’, analog=False, output=’ba’)

输入参数：
N:滤波器的阶数
Wn：归一化截止频率。计算公式Wn=2*截止频率/采样频率。（注意：根据采样定理，采样频率要大于两倍的信号本身最大的频率，才能还原信号。截止频率一定小于信号本身最大的频率，所以Wn一定在0和1之间）。当构造带通滤波器或者带阻滤波器时，Wn为长度为2的列表。
btype : 滤波器类型{‘lowpass’, ‘highpass’, ‘bandpass’, ‘bandstop’},
output : 输出类型{‘ba’, ‘zpk’, ‘sos’},

输出参数：
b，a: IIR滤波器的分子（b）和分母（a）多项式系数向量。output=’ba’
z,p,k: IIR滤波器传递函数的零点、极点和系统增益. output= ‘zpk’
sos: IIR滤波器的二阶截面表示。output= ‘sos’

低通滤波

scipy.signal.butter

butterworth 滤波

import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import scipy
from scipy import signal, fft
from scipy.signal import cheby1

# 显示中文标签
plt.rcParams['font.sans-serif']=['SimHei']
# 显示负号
plt.rcParams['axes.unicode_minus']=False
#有中文使用 u'中文内容'

def test():
    order = 5
    sampling_freq = 30
    cutoff_freq = 0.707
    sampling_duration = 5
    number_of_samples = sampling_freq * sampling_duration
    time = np.linspace(0, sampling_duration, number_of_samples, endpoint=False)
    signal = np.sin(2 * np.pi * time) + 0.5 * np.cos(6 * 2 * np.pi * time) + 1.5 * np.sin(9 * 2 * np.pi * time)
    normalized_cutoff_freq = 2 * cutoff_freq / sampling_freq
    numerator_coeffs, denominator_coeffs = scipy.signal.butter(order, normalized_cutoff_freq)
    filtered_signal = scipy.signal.lfilter(numerator_coeffs, denominator_coeffs, signal)
    plt.plot(time, signal, 'b-', label='signal')
    plt.plot(time, filtered_signal, 'g-', linewidth=2, label='filtered signal')
    plt.legend()
    plt.show()

if __name__ == '__main__':
    test()

这里假设采样频率为1000hz,信号本身最大的频率为500hz，要滤除400hz以上频率成分，即截至频率为400hz,则wn=2*400/1000=0.8。Wn=0.8

from scipy import signal
 
b, a = signal.butter(8, 0.8, 'lowpass')   #配置滤波器 8 表示滤波器的阶数
filtedData = signal.filtfilt(b, a, data)  #data为要过滤的信号

Chebyshev filter

## 切比雪夫1型低通滤波器，增益为1.0, 截止频率为 0.4
B, A = signal.cheby1(5, 1.0, 0.4)


fs = 1000;   % Hz采样频率
wp = 55/(fs/2);  %通带截止频率,取50~100中间的值,并对其归一化
ws = 90/(fs/2);  %阻带截止频率,取50~100中间的值,并对其归一化
Rp = 3;           %通带允许最大衰减为 db
Rs = 40;          %阻带允许最小衰减为 db
[n,Wn]=cheb1ord(Wp,Ws,Rp,Rs); % 获取阶数和截止频率
[b,a]=cheby1(n,Rp,Wn, 'low');  %获得转移函数系数
dataOut = filter(b,a,dataIn);          %信号滤波运算

高通滤波

这里假设采样频率为1000hz,信号本身最大的频率为500hz，要滤除100hz以下频率成分，即截至频率为100hz,则wn=2*100/1000=0.2。Wn=0.2

from scipy import signal
 
b, a = signal.butter(8, 0.2, 'highpass')   #配置滤波器 8 表示滤波器的阶数
filtedData = signal.filtfilt(b, a, data)  #data为要过滤的信号

带通滤波

这里假设采样频率为1000hz,信号本身最大的频率为500hz，要滤除100hz以下，400hz以上频率成分，即截至频率为100，400hz,则wn1=2100/1000=0.2，Wn1=0.2； wn2=2400/1000=0.8，Wn2=0.8。Wn=[0.02,0.8]

from scipy import signal
 
b, a = signal.butter(8, [0.2,0.8], 'bandpass')   #配置滤波器 8 表示滤波器的阶数
filtedData = signal.filtfilt(b, a, data)  #data为要过滤的信号

带阻滤波

这里假设采样频率为1000hz,信号本身最大的频率为500hz，要滤除100hz以上，400hz以下频率成分，即截至频率为100，400hz,则wn1=2100/1000=0.2，Wn1=0.2； wn2=2400/1000=0.8，Wn2=0.8。Wn=[0.02,0.8]，和带通相似，但是带通是保留中间，而带阻是去除。

from scipy import signal
 
b, a = signal.butter(8, [0.2,0.8], 'bandstop')   #配置滤波器 8 表示滤波器的阶数
filtedData = signal.filtfilt(b, a, data)  #data为要过滤的信号

Applying filter in scipy.signal: Use lfilter or filtfilt?

stackexchange

filtfilt is zero-phase filtering, which doesn’t shift the signal as it filters. Since the phase is zero at all frequencies, it is also linear-phase. Filtering backwards in time requires you to predict the future, so it can’t be used in “online” real-life applications, only for offline processing of recordings of signals.
lfilter is causal forward-in-time filtering only, similar to a real-life electronic filter. It can’t be zero-phase. It can be linear-phase (symmetrical FIR), but usually isn’t. Usually it adds different amounts of delay at different frequencies.

from __future__ import division, print_function
import numpy as np
from numpy.random import randn
from numpy.fft import rfft
from scipy import signal
import matplotlib.pyplot as plt

b, a = signal.butter(4, 0.03, analog=False)

# Show that frequency response is the same
impulse = np.zeros(1000)
impulse[500] = 1

# Applies filter forward and backward in time
imp_ff = signal.filtfilt(b, a, impulse)

# Applies filter forward in time twice (for same frequency response)
imp_lf = signal.lfilter(b, a, signal.lfilter(b, a, impulse))

plt.subplot(2, 2, 1)
plt.semilogx(20*np.log10(np.abs(rfft(imp_lf))))
plt.ylim(-100, 20)
plt.grid(True, which='both')
plt.title('lfilter')

plt.subplot(2, 2, 2)
plt.semilogx(20*np.log10(np.abs(rfft(imp_ff))))
plt.ylim(-100, 20)
plt.grid(True, which='both')
plt.title('filtfilt')

sig = np.cumsum(randn(800))  # Brownian noise
sig_ff = signal.filtfilt(b, a, sig)
sig_lf = signal.lfilter(b, a, signal.lfilter(b, a, sig))
plt.subplot(2, 1, 2)
plt.plot(sig, color='silver', label='Original')
plt.plot(sig_ff, color='#3465a4', label='filtfilt')
plt.plot(sig_lf, color='#cc0000', label='lfilter')
plt.grid(True, which='both')
plt.legend(loc="best")
plt.show()