Python 基于FIR實現(xiàn)Hilbert濾波器求信號包絡詳解
在通信領域,可以通過希爾伯特變換求解解析信號,進而求解窄帶信號的包絡。
實現(xiàn)希爾伯特變換有兩種方法,一種是對信號做FFT,單后只保留單邊頻譜,在做IFFT,我們稱之為頻域方法;另一種是基于FIR根據(jù)傳遞函數(shù)設計一個希爾伯特濾波器,我們稱之為時域方法。
# -*- coding:utf8 -*-# @TIME : 2019/4/11 18:30# @Author : SuHao# @File : hilberfilter.pyimport scipy.signal as signalimport numpy as npimport librosa as libimport matplotlib.pyplot as pltimport time# from preprocess_filter import *# 讀取音頻文件ex = ’....數(shù)據(jù)集2pre2012bfluteBassFlute.ff.C5B5.aiff’time_series, fs = lib.load(ex, sr=None, mono=True, res_type=’kaiser_best’)# 生成一個chirp信號# duration = 2.0# fs = 400.0# samples = int(fs*duration)# t = np.arange(samples) / fs# time_series = signal.chirp(t, 20.0, t[-1], 100.0)# time_series *= (1.0 + 0.5 * np.sin(2.0*np.pi*3.0*t) )def hilbert_filter(x, fs, order=201, pic=None): ’’’ :param x: 輸入信號 :param fs: 信號采樣頻率 :param order: 希爾伯特濾波器階數(shù) :param pic: 是否繪圖,bool :return: 包絡信號 ’’’ co = [2*np.sin(np.pi*n/2)**2/np.pi/n for n in range(1, order+1)] co1 = [2*np.sin(np.pi*n/2)**2/np.pi/n for n in range(-order, 0)] co = co1+[0]+ co # out = signal.filtfilt(b=co, a=1, x=x, padlen=int((order-1)/2)) out = signal.convolve(x, co, mode=’same’, method=’direct’) envolope = np.sqrt(out**2 + x**2) if pic is not None: w, h = signal.freqz(b=co, a=1, worN=2048, whole=False, plot=None, fs=2*np.pi) fig, ax1 = plt.subplots() ax1.set_title(’hilbert filter frequency response’) ax1.plot(w, 20 * np.log10(abs(h)), ’b’) ax1.set_ylabel(’Amplitude [dB]’, color=’b’) ax1.set_xlabel(’Frequency [rad/sample]’) ax2 = ax1.twinx() angles = np.unwrap(np.angle(h)) ax2.plot(w, angles, ’g’) ax2.set_ylabel(’Angle (radians)’, color=’g’) ax2.grid() ax2.axis(’tight’) # plt.savefig(pic + ’hilbert_filter.jpg’) plt.show() # plt.clf() # plt.close() return envolopestart = time.time()env0 = hilbert_filter(time_series, fs, 81, pic=True)end = time.time()a = end-startprint(a)plt.figure()ax1 = plt.subplot(211)plt.plot(time_series)ax2 = plt.subplot(212)plt.plot(env0)plt.xlabel(’time’)plt.ylabel(’mag’)plt.title(’envolope of music by FIR n time:%.3f’%a)plt.tight_layout()start = time.time()# 使用scipy庫函數(shù)實現(xiàn)希爾伯特變換env = np.abs(signal.hilbert(time_series))end = time.time()a = end-startprint(a)plt.figure()ax1 = plt.subplot(211)plt.plot(time_series)ax2 = plt.subplot(212)plt.plot(env)plt.xlabel(’time’)plt.ylabel(’mag’)plt.title(’envolope of music by scipy n time:%.3f’%a)plt.tight_layout()plt.show()
使用chirp信號對兩種方法進行比較
FIR濾波器的頻率響應
使用音頻信號對兩種方法進行比較
由于音頻信號時間較長,采樣率較高,因此離散信號序列很長。使用頻域方法做FFT和IFFT要耗費比較長的時間;然而使用時域方法只是和濾波器沖擊響應做卷積,因此運算速度比較快。結果對比如下:
頻域方法結果
時域方法結果
由此看出,時域方法耗費時間要遠小于頻域方法。
以上這篇Python 基于FIR實現(xiàn)Hilbert濾波器求信號包絡詳解就是小編分享給大家的全部內(nèi)容了,希望能給大家一個參考,也希望大家多多支持好吧啦網(wǎng)。
相關文章:
