假设有一段10kHz的语言，现需要对2~3kHz之间的语言信号进行提取，要求1.5kHz及3.5kHz以上的频率需要有40dB的衰减
1、求数字频率指标
通带下边频：
w
p
l
=
2
∗
π
∗
f
p
l
/
f
s
=
0.4
π
w_{pl}=2*\pi *f_{pl}/f_s=0.4\pi
wpl=2∗π∗fpl/fs=0.4π
通带上边频：
w
p
h
=
2
∗
π
∗
f
p
h
/
f
s
=
0.6
π
w_{ph}=2*\pi *f_{ph}/f_s=0.6\pi
wph=2∗π∗fph/fs=0.6π
下阻带上变频：
w
s
l
=
2
∗
π
∗
f
s
l
/
f
s
=
0.3
π
w_{sl}=2*\pi *f_{sl}/f_s=0.3\pi
wsl=2∗π∗fsl/fs=0.3π
上阻带下变频：
w
s
h
=
2
∗
π
∗
f
s
h
/
f
s
=
0.7
π
w_{sh}=2*\pi *f_{sh}/f_s=0.7\pi
wsh=2∗π∗fsh/fs=0.7π
2、选取窗函数
根据阻带衰减查表，可选汉宁窗,过度带宽
Δ
w
=
w
p
l
−
w
s
l
=
0.1
π
\Delta_w=w_{pl}-w_{sl}=0.1\pi
Δw=wpl−wsl=0.1π
由汉宁窗过度带宽确定阶数N
N
=
6.2
π
/
Δ
w
=
62
N=6.2\pi/\Delta_w=62
N=6.2π/Δw=62
取N为奇数N=63
a
=
(
N
−
1
)
/
2
a = (N-1)/2
a=(N−1)/2
因此窗函数：
w
(
n
)
=
1
2
[
1
−
c
o
s
(
2
π
n
a
)
]
w(n)=\frac{1}{2}[1-cos(\frac{2\pi n}{a})]
w(n)=21[1−cos(a2πn)]
3、求理想带通滤波器的单位脉冲响应
理想带通滤波器的截止频率：
w
c
l
=
(
w
p
l
+
w
s
l
)
/
2
w_{cl}=(w_{pl}+w_{sl})/2
wcl=(wpl+wsl)/2
w
c
h
=
(
w
p
h
+
w
s
h
)
/
2
w_{ch}=(w_{ph}+w_{sh})/2
wch=(wph+wsh)/2
理想带通滤波器的单位脉冲响应：
h d ( n ) = s i n [ w c h ∗ ( n − a ) ] − s i n [ w c l ∗ ( n − a ) ] π ∗ ( n − a ) h_d(n)=\frac{sin[w_{ch}*(n-a)]-sin[w_{cl}*(n-a)]}{\pi*(n-a)} hd(n)=π∗(n−a)sin[wch∗(n−a)]−sin[wcl∗(n−a)]
4、求FIR滤波参数
h ( n ) = h d ( n ) w ( n ) h(n)=h_d(n)w(n) h(n)=hd(n)w(n)
5、算法仿真
import numpy as np
import matplotlib.pyplot as plt
from scipy import signal
from scipy.fftpack import fft,ifft
#采样为10Khz
#1.5khz以下及3.5khz以上至少40db的衰减
f_pl = 2000
f_sl = 1500
f_sh = 3500
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,a-1,a)
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
h_d = (np.sin(W_ch*(n-a))-np.sin(W_cl*(n-a)))/(np.pi*(n-a))
#h_d = (np.sin(0.65*(n-31)*np.pi)-np.sin(0.35*(n-31)*np.pi))/(np.pi*(n-31))
h_c = h_d*w_n
numtaps=63
array = [0]*63
for i in range(31):
array[i] = h_d[i]
array[60-i] = -h_d[i]
plt.figure(1)
#plt.plot(h_d,'b')
plt.title('滤波系数') # 定义标题
plt.plot(array,'g')
#plt.plot(n,h_c,'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(32,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,'r') #显示原始信号的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,'r')
plt.subplot(224)
plt.title('滤波后的信号FFT') # 定义标题
plt.plot(xf,yf1_1,'r') #显示原始信号的FFT模值