假设有一段10kHz的语言，现需要对2~3kHz之间的语言信号进行提取，要求1.5kHz及3.5kHz以上的频率需要有40dB的衰减

1、求数字频率指标

通带下边频：
$$w_{pl}=2\ast \pi \ast f_{pl}/f_s=0.4\pi$$
通带上边频：
$$w_{ph}=2\ast \pi \ast f_{ph}/f_s=0.6\pi$$
下阻带上变频：
$$w_{sl}=2\ast \pi \ast f_{sl}/f_s=0.3\pi$$
上阻带下变频：
$$w_{sh}=2\ast \pi \ast f_{sh}/f_s=0.7\pi$$

2、选取窗函数

根据阻带衰减查表，可选汉宁窗,过度带宽$${\Delta }_w=w_{pl}-w_{sl}=0.1\pi$$
由汉宁窗过度带宽确定阶数N
$$N=6.2\pi /{\Delta }_w=62$$
取N为奇数N=63
$$a=(N-1)/2$$
因此窗函数：
$$w(n)=\frac{1}{2}[1-cos(\frac{2\pi n}{a})]$$

3、求理想带通滤波器的单位脉冲响应

理想带通滤波器的截止频率：
$$w_{cl}=(w_{pl}+w_{sl})/2$$
$$w_{ch}=(w_{ph}+w_{sh})/2$$
理想带通滤波器的单位脉冲响应：

$$h_d(n)=\frac{sin[w_{ch}\ast (n-a)]-sin[w_{cl}\ast (n-a)]}{\pi \ast (n-a)}$$

4、求FIR滤波参数

$$h(n)=h_d(n)w(n)$$

5、算法仿真

import numpy as np
import matplotlib.pyplot as plt
from scipy import signal
from scipy.fftpack import fft,ifft
from decimal import Decimal
import matplotlib
matplotlib.rcParams['font.sans-serif'] = ['SimHei']
matplotlib.rcParams['font.family']='sans-serif'

 
class filter:
    def __init__(self,h):
        self.order=len(h)
        self.h=h
        self.output=[]
    def FIR_Filter(self,vi):
        for i in range(len(vi)):
            sum=0
            if i < self.order:
                for j in range(i):
                    sum=sum + self.h[j]*vi[i-j]
            else:      
                for j in range(self.order):
                    sum=sum + self.h[j]*vi[i-j]
                
            self.output.append(sum)   
        return self.output
        
       

#采样为10Khz
#1.5khz以下及3.5khz以上至少40db的衰减

f_sl = 1500
f_sh = 3500
f_pl = 2000
f_ph = 3000
f_s  = 10000
#通带下边频
W_pl = 2*np.pi*f_pl/f_s

W_ph = 2*np.pi*f_ph/f_s

W_sl = 2*np.pi*f_sl/f_s
W_sh = 2*np.pi*f_sh/f_s

W_D = W_pl - W_sl

N = 6.2*np.pi/(W_D)
if N%2==0:
    N=N+1
print(N)
a = (N-1)/2
n=np.linspace(0,N-1,N)

R_n =  1
#汉宁窗口函数
w_n = 0.5*(1-np.cos(2*np.pi*n/(N-1)))

W_cl = (W_pl+W_sl)/2
W_ch = (W_ph+W_sh)/2
#用一个靠近a的小数将a值替换掉，避免出现除0的情况
a=30.9999999999
h_d  = (np.sin(W_ch*(n-a))-np.sin(W_cl*(n-a)))/(2*np.pi*(n-a))
h_c = h_d*w_n

numtaps=63

array= h_c
plt.figure(1)
yy_1=fft(array)                     #快速傅里叶变换
yf_1=abs(fft(array))                # 取模
yf1_1=abs(fft(array))/((len(array)/2))           #归一化处理
yf2_1 = yf1_1[range(int(len(array)/2))]  #由于对称性，只取一半区间
#plt.plot(h_d,'b')
plt.subplot(221)
plt.title('滤波系数')  # 定义标题
plt.plot(array,'g')
plt.plot(h_c,'K')
plt.subplot(222)
plt.title('滤波系数FFT')  # 定义标题
plt.plot(yf2_1,'r')
plt.show()

x=np.linspace(0,1,f_s)
signal_array = np.sin(2*np.pi*2000*x)
for i in range(10):
    if 1000*i != 2000:
        signal_array+=np.sin(2*np.pi*1000*x*i)#+np.sin(2*np.pi*175*x)+np.sin(2*np.pi*350*x)+np.sin(2*np.pi*500*x)
plt.figure(2)        
Weight = array
FIR_filter=filter(Weight)
output = FIR_filter.FIR_Filter(signal_array)   
y=  signal_array
xf = np.arange(len(y)) 
yy=fft(y)                     #快速傅里叶变换
yf=abs(fft(y))                # 取模
yf1=abs(fft(y))/((len(x)/2))           #归一化处理
yf2 = yf1[range(int(len(x)/2))]  #由于对称性，只取一半区间
plt.subplot(221)
plt.title('原始信号')  # 定义标题
plt.plot(xf,signal_array,'b') #显示原始信号的FFT模值

plt.subplot(222)
plt.title('原始信号FFT')  # 定义标题
plt.plot(xf,yf1,'r') #显示原始信号的FFT模值

yy_1=fft(output)                     #快速傅里叶变换
yf_1=abs(fft(output))                # 取模
yf1_1=abs(fft(output))/((len(x)/2))           #归一化处理
yf2_1 = yf1_1[range(int(len(x)/2))]  #由于对称性，只取一半区间
plt.subplot(223)
plt.title('滤波后的信号')  # 定义标题
plt.plot(xf,output,'b') 
plt.subplot(224)
plt.title('滤波后的信号FFT')  # 定义标题
plt.plot(xf,yf1_1,'r') #显示原始信号的FFT模值

        

6、算法结果