利用Python实现IIR滤波器

1.1 直接I型

1.2 算法实现

class filter:
    def __init__(self,):
        pass
   

    def IIR_Filter_I(self,input_array,a_weight,b_weight,Scale_Factors):
        output_array=[0,0,0]
        y_0 = 0
        y_1 = 0
        y_2 = 0
        x_0 = input_array[0]
        x_1 = 0 
        x_2 = 0
        for index_vi in range(1,len(input_array)):
            x_0 = x_1
            x_1 = x_2
            x_2 = input_array[index_vi]/Scale_Factors
            y_0 = y_1
            y_1 = y_2
            y_2 = (x_0 * b_weight[0] + x_1 *b_weight[1] + x_2*b_weight[2] + y_0 * a_weight[1] + y_1 *a_weight[2])

            output_array.append(y_2)  
        return output_array
    
    
    def IIR_Filter_II(self,vi):
        pass
    

1.3 数据仿真

以下代码实现一个二阶高通滤波器，蓝色的是原数据，红色的是滤波后的数据：

对原数据和滤波后的数据进行FFT分析：

可见低于200Hz的信号被滤掉，红线是原数据的频谱图，蓝线对应滤波后的频谱图。
采用的是fdatool工具，生成的是截止频率是200Hz的滤波参数。
完整代码如下：

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

a_weight=[1 , -0.369527377351242,  0.195815712655833]
b_weight=[1 , -2 , 1]
Scale_Factors =0.391335772501769  
class filter:
    def __init__(self,):
        pass
   

    def IIR_Filter_I(self,input_array,a_weight,b_weight,Scale_Factors):
        output_array=[0,0,0]
        y_0 = 0
        y_1 = 0
        y_2 = 0
        x_0 = input_array[0]
        x_1 = 0 
        x_2 = 0
        for index_vi in range(1,len(input_array)):
            
            y = (x_0 * b_weight[0] + x_1 *b_weight[1] + x_2*b_weight[2] + y_0 * a_weight[0] + y_1 *a_weight[1]+ y_2 *a_weight[2])*Scale_Factors
            x_2=x_1
            x_1=x_0
            x_0=input_array[index_vi]
            
            y_2=y_1
            y_1=y_0
            y_0=y
           
            output_array.append(y*Scale_Factors)  
        return output_array
    
    
    def IIR_Filter_II(self,vi):
        pass
    
    
#采样点选择1400个，因为设置的信号频率分量最高为600Hz，根据采样定理知采样频率要大于信号频率2倍，所以这里设置采样频率为1400Hz（即一秒内有1400个采样点）
x=np.linspace(0,1,1200)
#设置需要采样的信号，频率分量有180，390和600
signal=np.sin(2*np.pi*50*x) + np.sin(2*np.pi*150*x) +  np.sin(2*np.pi*400*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（双边频率范围）

IIR_Output =  IIR_filter.IIR_Filter_I(signal,a_weight,b_weight,Scale_Factors)
IIR_Output_Size = len(IIR_Output)
plt.figure(1)
plt.plot(signal[0:50],'b') #显示原始信号的FFT模值
plt.plot(IIR_Output[0:50],'r') #显示原始信号的FFT模值

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:1200],'r') #显示原始信号的FFT模值
plt.plot(IIR_Output_fft_abs_size,IIR_Output_fft_abs_norm,'g') #显示原始信号的FFT模值
plt.show()

2.1直接II型

2.2 算法实现

class filter:
    def __init__(self,):
        pass
   

    def IIR_Filter_II(self,input_array,a_weight,b_weight,Scale_Factors):
        output_array=[0]
  
        x_0 = input_array[0]

        
        w_0 = 0
        w_1 = 0
        w_2 = 0
        
        for index_vi in range(0,len(input_array)):
            x_0 = input_array[index_vi]
            w_n = x_0 - ( a_weight[0]*w_0 + a_weight[1]*w_1 + a_weight[2]*w_2)
          
            y = b_weight[0]*w_0 + b_weight[1]*w_1 + b_weight[1]*w_2 
            #w_0 = w_1
            w_2 = w_1
            w_1 = w_0
            w_0 = w_n
            output_array.append(y*Scale_Factors )  
        return output_array
    

2.3 数据仿真

略