Python 基于FIR实现Hilbert滤波器求信号包络详解

作者:qq7835144@163.com 时间:2023-07-13 01:31:47 

在通信领域,可以通过希尔伯特变换求解解析信号,进而求解窄带信号的包络。

实现希尔伯特变换有两种方法,一种是对信号做FFT,单后只保留单边频谱,在做IFFT,我们称之为频域方法;另一种是基于FIR根据传递函数设计一个希尔伯特滤波器,我们称之为时域方法。


# -*- coding:utf8 -*-
# @TIME   : 2019/4/11 18:30
# @Author  : SuHao
# @File   : hilberfilter.py

import scipy.signal as signal
import numpy as np
import librosa as lib
import matplotlib.pyplot as plt
import time
# from preprocess_filter import *

# 读取音频文件
ex = '..\\..\\数据集2\\pre2012\\bflute\\BassFlute.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: 希尔伯特滤波器阶数
 :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 envolope

start = time.time()
env0 = hilbert_filter(time_series, fs, 81, pic=True)
end = time.time()
a = end-start
print(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库函数实现希尔伯特变换
env = np.abs(signal.hilbert(time_series))
end = time.time()
a = end-start
print(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滤波器的频率响应

Python 基于FIR实现Hilbert滤波器求信号包络详解

使用音频信号对两种方法进行比较

由于音频信号时间较长,采样率较高,因此离散信号序列很长。使用频域方法做FFT和IFFT要耗费比较长的时间;然而使用时域方法只是和滤波器冲击响应做卷积,因此运算速度比较快。结果对比如下:

频域方法结果

Python 基于FIR实现Hilbert滤波器求信号包络详解

时域方法结果

Python 基于FIR实现Hilbert滤波器求信号包络详解

由此看出,时域方法耗费时间要远小于频域方法。

来源:https://blog.csdn.net/qq7835144/article/details/102811793

标签:Python,滤波器,信号,包络
0
投稿

猜你喜欢

  • 一段压缩MS SQLServer日志的语句

    2024-01-14 23:42:07
  • php实现文件下载更能介绍

    2023-08-18 14:30:21
  • FrontPage 2002应用技巧四则

    2008-08-17 10:57:00
  • 搞清楚 Python traceback的具体使用方法

    2022-10-04 07:00:15
  • python中实现词云图的示例

    2021-08-04 22:11:33
  • 基于Python实现图片一键切割九宫格的工具

    2022-07-05 00:22:31
  • SQL查询之字段是逗号分隔开的数组如何查询匹配数据问题

    2024-01-21 22:32:08
  • 打开电脑上的QQ的python代码

    2022-08-18 04:21:28
  • Python rstrip()方法实例详解

    2021-05-16 10:05:40
  • python数据处理 根据颜色对图片进行分类的方法

    2022-02-27 08:41:37
  • python DataFrame 修改列的顺序实例

    2023-07-21 12:55:40
  • Django显示可视化图表的实践

    2023-04-13 02:42:38
  • 对python 树状嵌套结构的实现思路详解

    2022-02-04 15:45:06
  • python linecache读取行更新的实现

    2021-01-26 01:33:06
  • python实现的AES双向对称加密解密与用法分析

    2022-08-11 00:42:32
  • PHP使用flock实现文件加锁的方法

    2023-10-29 21:26:59
  • Node.js查找当前目录下文件夹实例代码

    2024-04-16 09:50:32
  • Mysql 数据库死锁过程分析(select for update)

    2024-01-23 02:57:26
  • asp如何正确理解和使用Command、Connection和 Recordset三个对象?

    2010-06-28 18:23:00
  • asp中把数据表映射成ajax可调用的json格式的方法

    2010-01-22 15:27:00
  • asp之家 网络编程 m.aspxhome.com