大家好,又見面了,我是你們的朋友全棧君。
假設有一段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
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、算法結果
釋出者:全棧程式員棧長,轉載請注明出處:https://javaforall.cn/138893.html原文連結:https://javaforall.cn