基于python实现IIR带通滤波器（该文存在问题，仅供读者参考，不能当做学习范例使用）

该文存在问题，仅供读者参考，不能当做学习范例使用，感兴趣的读者可访问：
Python 实现巴特沃斯滤波器
巴特沃斯滤波器是计算IIR参数的一种方法

1、一阶IIR算法

def IIR_I(self,input_array,a_weight,b_weight,Scale_Factors):
        output_array = []
        for index_vi in range(0,len(input_array)):
            xv[0] = xv[1]
            xv[1] = xv[2]
            xv[2] = input_array[index_vi]/Scale_Factors
            
            yv[0] = yv[1]
            yv[1] = yv[2]
            yv[2] = xv[0] * b_weight[0] + xv[1]* b_weight[1]  + xv[2]* b_weight[2] +  a_weight[1] * yv[0] +  a_weight[2] * yv[1]
            yv[2] = yv[2]*Scale_Factors
            output_array.append(yv[2])  
        return output_array
    

2、利用fdatool生成滤波参数

 1.0f,  0.0f,  -1.0f,    -1.1276518720541668f,   -0.47001314508753411f,       
  1.0f,  0.0f,  -1.0f,    0.77495305804604886f,  -0.36707750055668387f   

3、算法实现

# -*- coding: utf-8 -*-
"""
Created on Wed Mar 17 14:04:40 2021

@author: zhanghui
"""

import numpy as np
import matplotlib.pyplot as plt
from scipy.fftpack import fft,ifft


NZEROS = 2
NPOLES = 2
xv=[0] * (NZEROS+1)
yv=[0] * (NPOLES+1)
class filter:
    def __init__(self,):
        pass

    def IIR_I(self,input_array,a_weight,b_weight,Scale_Factors):
        
       
        output_array = []
        for index_vi in range(0,len(input_array)):
            xv[0] = xv[1]
            xv[1] = xv[2]
            xv[2] = input_array[index_vi]/Scale_Factors
            
            yv[0] = yv[1]
            yv[1] = yv[2]
            yv[2] = xv[0] * b_weight[0] + xv[1]* b_weight[1]  + xv[2]* b_weight[2] +  a_weight[1] * yv[0] +  a_weight[2] * yv[1]
            yv[2] = yv[2]*Scale_Factors
            output_array.append(yv[2])  
        return output_array
    
  
#采样点选择1400个，因为设置的信号频率分量最高为600Hz，根据采样定理知采样频率要大于信号频率2倍，所以这里设置采样频率为1400Hz（即一秒内有1400个采样点）
fs = 1400
x=np.linspace(0,1,fs)
#设置需要采样的信号，频率分量有180，390和600
signal=np.sin(2*np.pi*50*x)+np.sin(2*np.pi*300*x)# + np.sin(2*np.pi*500*x)#+np.sin(2*np.pi*400*x) + np.sin(2*np.pi*500*x) +  np.sin(2*np.pi*600*x) #+  np.sin(2*np.pi*6500*x) +  np.sin(2*np.pi*7000*x)

signal_fft=fft(signal)                     #快速傅里叶变换
signal_fft_abs=abs(fft(signal))                # 取模
signal_fft_abs_norm=abs(fft(signal))/((len(signal)/2))           #归一化处理
signal_fft_abs_norm_half= signal_fft_abs_norm[range(int(len(signal)/2))]  #由于对称性，只取一半区间
signal_fft_abs_size =np.arange(len(signal))        # 频率
IIR_filter=filter()
#混合波的FFT（双边频率范围）


a_weight=[1,   0.77495305804604886,  -0.36707750055668387]
b_weight=[1 , 0 , -1]
Scale_Factors = 0.55815658576077365

IIR_Output =  IIR_filter.IIR_I(signal,a_weight,b_weight,Scale_Factors)

a_weight=[1 ,   -1.1276518720541668,   -0.47001314508753411]
b_weight=[1 , 0 , -1]
Scale_Factors = 0.55815658576077365 

IIR_Output =  IIR_filter.IIR_I(IIR_Output,a_weight,b_weight,Scale_Factors)
#a_weight=[1 , -1.996481382567631,  0.996532750562765]
#IIR_Output =  IIR_filter.filterloop(IIR_Output,a_weight)
#IIR_Output =  IIR_filter.filterloop(signal)

IIR_Output_Size = len(IIR_Output)
plt.figure(1)
plt.plot(signal[0:500],'g') 
plt.plot(IIR_Output[0:500],'b') 

IIR_Output = np.array(IIR_Output)
IIR_Output_fft=fft(IIR_Output)                     #快速傅里叶变换
IIR_Output_fft_abs=abs(fft(IIR_Output))                # 取模
IIR_Output_fft_abs_norm=abs(fft(IIR_Output))/((len(IIR_Output)/2))           #归一化处理
IIR_Output_fft_abs_half = IIR_Output_fft_abs_norm[range(int(len(IIR_Output)/2))]  #由于对称性，只取一半区间
IIR_Output_fft_abs_size = np.arange(len(IIR_Output_fft_abs_norm))        # 频率

plt.figure(2)
plt.plot(signal_fft_abs_size,signal_fft_abs_norm[0:fs],'r') #显示原始信号的FFT模值
plt.plot(IIR_Output_fft_abs_size,IIR_Output_fft_abs_norm,'g') 
plt.show()

4、实验结果