基于python3实现的数字滤波器[0]

Posted on by

本篇是基于部署jupyterhub后实现的python工程计算脚本。

  • 低通滤波器
import numpy as np
from matplotlib import pyplot as plt
f=1 #信号频率Hz
t=np.linspace(0,5*1/f,5*int(1/f*1000)) #时间
x=np.sin(2*np.pi*f*t)  #生成信号曲线


#滤波算法
a0=0.9691
b0=0.01547
b1=0.01547
n=len(t)
y=np.zeros((1,n),dtype=float)

#截至频率5Hz，采样频率1KHz
for i in range(0,n-1):
    y[:,i+1]=b1*x[i+1]+b0*x[i]+a0*y[:,i]

#绘制滤波后信号
plt.plot(t,x,'r')
plt.plot(t,y[0,:],'--')
plt.legend(['input','output'])
plt.grid()
plt.show()

1Hz效果

6Hz效果如下：

  • 高通滤波器
import numpy as np
from matplotlib import pyplot as plt
f=10 #信号频率Hz
t=np.linspace(0,5*1/f,5*int(1/f*1000)) #时间
n=len(t)
x=np.sin(2*np.pi*f*t)  #生成信号曲线


#滤波算法
a0=0.9691
b0=-0.9845
b1=0.9845
y=np.zeros((1,n),dtype=float)

#截至频率5Hz，采样频率1KHz
for i in range(0,n-1):
    y[:,i+1]=b1*x[i+1]+b0*x[i]+a0*y[:,i]

#绘制滤波后信号
plt.plot(t,x,'r')
plt.plot(t,y[0,:],'--')
plt.legend(['input','output'])
plt.grid()
plt.show()

10Hz

1Hz效果如下：

  • 带通滤波器
import numpy as np
from matplotlib import pyplot as plt
f=1 #信号频率Hz，1 10 60
t=np.linspace(0,5*1/f,5*int(1/f*1000)) #时间
n=len(t)
x=np.sin(2*np.pi*f*t)  #生成信号曲线


#滤波算法
a0=-0.7164
a1=1.712
b0=-0.1346
b1=0
b2=0.1346
y=np.zeros((1,n),dtype=float)

#频率2.6-50Hz带通滤波器，采样频率1KHz
for i in range(0,n-2):
    y[:,i+2]=b2*x[i+2]+b0*x[i]+a1*y[:,i+1]+a0*y[:,i]

#绘制滤波后信号
plt.plot(t,x,'r')
plt.plot(t,y[0,:],'--')
plt.legend(['input','output'])
plt.grid()
plt.show()

1Hz效果如下：

10Hz效果

60Hz效果

  • butterworth滤波器
import numpy as np
from matplotlib import pyplot as plt
f=1 #信号频率Hz ,1,20
t=np.linspace(0,5*1/f,5*int(1/f*1000)) #时间
n=len(t)
x=10*np.sin(2*np.pi*f*t)  #生成信号曲线

#滤波算法
a0=-0.84860184
a1=3.5336803
a2=-5.5209723
a3=3.8358795
b0=0.89732853/pow(10,6)
b1=3.5893141/pow(10,6)
b2=5.3839712/pow(10,6)
b3=-3.58931411/pow(10,6)
b4=0.89732853/pow(10,6)

y=np.zeros((1,n),dtype=float)

#10hz频率进行滤波，采样频率1KHz
#为啥需要2倍？
for i in range(0,n-4):
    y[:,i+4]=2*(b4*x[i+4]+b3*x[i+3]+b2*x[i+2]+b1*x[i+1]+b0*x[i]) \
              +a3*y[:,i+3]+a2*y[:,i+2]+a1*y[:,i+1]+a0*y[:,i]


#绘制滤波后信号
plt.plot(t,x,'r')
plt.plot(t,y[0,:],'--')
plt.legend(['input','output'])
plt.grid()
plt.show()

1Hz效果：

20Hz效果

  • 贝塞尔滤波器
import numpy as np
from matplotlib import pyplot as plt
f=1 #信号频率Hz ,1,10
t=np.linspace(0,10*1/f,10*int(1/f*1000)) #时间
n=len(t)
x=10*np.sin(2*np.pi*f*t)  #生成信号曲线

#滤波算法
a0=0.9759648
a1=-2.9516865
a2=2.9757207
b0=0.12349019/pow(10,6)
b1=0.37047056/pow(10,6)
b2=0.37047056/pow(10,6)
b3=0.1234909/pow(10,6)

y=np.zeros((1,n),dtype=float)

#10hz频率进行滤波，采样频率1KHz
for i in range(0,n-3):
    y[:,i+3]=(b3*x[i+3]+b2*x[i+2]+b1*x[i+1]+b0*x[i]) \
              +a2*y[:,i+2]+a1*y[:,i+1]+a0*y[:,i]


#绘制滤波后信号
plt.plot(t,x,'r')
plt.plot(t,y[0,:],'--')
plt.legend(['input','output'])
plt.grid()
plt.show()

1Hz效果

20Hz效果