python - 在具有初始条件的连续块上过滤后重新采样（以避免不连续性）
问题描述
我实际上是在 python 上开发一个实时图形均衡器。我正在使用pyaudio模块scipy，，numpy。我的均衡器基于从 25 Hz 到 20 kHz（即 30 个频段）的第三倍频程滤波器组。该滤波器组将输入信号分成 30 个滤波信号（以每三个倍频程的中心频率为中心）。此外，流式传输是逐块实现的（使用 pyaudio 和回调方法）。
我filtfilt从scipy.signal模块中使用过，但每个块之间有一些不连续性（一些可听见的咔嗒声）。因此，在连续时间帧上应用 IIR 滤波器时，我遵循了连续性问题，它适用于高频。
但是对于低频，我需要按照以下步骤操作：
1）对输入信号进行下采样（以保持良好的清晰度过滤器）；
2）过滤lfilter_zi以保持每个块之间的连续性（用于流式传输）；
3) 对滤波后的信号进行上采样。
我的问题是上采样，因为这会破坏每个块之间的连续性（见下图） 下采样和上采样信号时 2 个块之间的不连续性（此处为 1000Hz 的正弦波）
另外，这是我的第三倍频程滤波器组代码：
# -*- coding: utf-8 -*-
"""
Created on Tue Feb 19 16:14:09 2019
@author: William
"""
from __future__ import division
import numpy as np
import scipy.signal as sc
def oct3_dsgn(fc, fs, type_filt='cheby2_bandpass', output='ba'):
"""
Calcul les coefficients B et A d'un filtre passe-bande pour chacune des
bandes de tiers d'octave entre 25 Hz et 20 kHz.
La fonction cheb2ord permet d'optimiser l'ordre et la bande passante du
filtre en fonction des fréquences de coupures et de leurs gains respectifs
(gpass, gstop) et de la fréquence d'échantillonnage.
"""
#------- Définition des fréquences inférieures et supérieures de la bande n
fc1 = fc / 2**(1/6)
fc2 = fc * 2**(1/6)
#------- Définition des fréquences centrales des bandes n-1 et n+1
fm1 = fc1 / 2**(1/6)
fm2 = fc2 * 2**(1/6)
#------- Définition des fréquences normalisées par rapport à f_nyquist
W1 = fc1/(fs/2)
W2 = fc2/(fs/2)
Wm1 = fm1/(fs/2)
Wm2 = fm2/(fs/2)
#------- Définition des filtres passe-bande, passe-bas et passe-haut
if type_filt == 'cheby2_bandpass':
gpass = 20*np.log10(np.sqrt(2)) # Équivalent à 3dB
gstop = 45
n_filt, Wn = sc.cheb2ord([W1, W2], [Wm1, Wm2], gpass, gstop)
if output=='ba':
B, A = sc.cheby2(n_filt, gstop, Wn, btype='band', output='ba')
elif output=='sos':
sos = sc.cheby2(n_filt, gstop, Wn, btype='band', output='sos')
elif type_filt == 'butter':
gpass = 20*np.log10(np.sqrt(2)) # Équivalent à 3dB
gstop = 30
n_filt, Wn = sc.buttord([W1, W2], [Wm1, Wm2], gpass, gstop)
if output=='ba':
B, A = sc.butter(n_filt, Wn, btype='band', output='ba')
elif output=='sos':
sos = sc.cheby2(n_filt, gstop, Wn, btype='band', output='sos')
elif type_filt == 'cheby2_low':
gpass = 20*np.log10(np.sqrt(2)) # Équivalent à 3dB
gstop = 45
n_filt, Wn = sc.cheb2ord(W2, Wm2, gpass, gstop)
if output=='ba':
B, A = sc.cheby2(n_filt, gstop, Wn, btype='low', output='ba')
elif output=='sos':
sos = sc.cheby2(n_filt, gstop, Wn, btype='low', output='sos')
elif type_filt == 'cheby2_high':
gpass = 20*np.log10(np.sqrt(2)) # Équivalent à 3dB
gstop = 45
n_filt, Wn = sc.cheb2ord(W1, Wm1, gpass, gstop)
if output=='ba':
B, A = sc.cheby2(n_filt, gstop, Wn, btype='high', output='ba')
elif output=='sos':
sos = sc.cheby2(n_filt, gstop, Wn, btype='high', output='sos')
if output == 'ba':
return B, A
elif output == 'sos':
return sos
def oct3_filter(signal, fs, gain_general = 0, gain_bande = np.zeros((30, )), plot=False):
"""
Calcul le signal filtré à partir du signal initial grâce une reconstruction
parfaite bande par bande. Le signal initial est filtré pour chacune des bandes.
Lorsque les fréquences sont trop basses, un sous-échantillonnage est opéré
pour gagner en résolution fréquentielle et en bande passante.
Pour la bande à 25 Hz, un filtre passe-bas est utilisé. La bande à 25 Hz
comprend donc toutes les fréquences inférieures ou égales à 25 Hz.
De même, pour la bande à 20 kHz un filtre passe-haut est utlisé. La bande à
20 kHz comprend donc toutes les fréquences au-dessus de 18 kHz.
"""
#------- Définition des fréquences centrales exactes en base 2
fc_oct3 = (1000) * (2**(1/3))**np.arange(-16, 14)
n_bande = len(fc_oct3)
n_signal = len(signal)
#------- Définition des matrices de stockage des signaux et des gains
signal_filt = np.zeros((n_signal, n_bande))
Gain = 10**(gain_bande/20)
#------- Affichage de la réponse des filtres 1/3 d'octaves
if plot == True:
import matplotlib.pyplot as plt
plt.figure(0, figsize=(10,5))
#------- Boucle sur les bandes de 10 Hz à 20 kHz
for ii in range(0, n_bande):
if ii == n_bande-1 :
## bande à 20 kHz
sos = oct3_dsgn(fc_oct3[ii], fs, type_filt='cheby2_high', output='sos')
signal_filt[:, ii] = Gain[ii] * sc.sosfilt(sos, signal)
## affichage de la réponse du filtre
if plot == True:
w, h = sc.freqz(sos)
plt.semilogx(w/2/np.pi*fs, 20*np.log10(abs(h)), 'k')
elif ii == 0 :
## bande à 25 Hz
n_decimate = 32
x = sc.decimate(signal, n_decimate)
sos = oct3_dsgn(fc_oct3[ii], fs//n_decimate, type_filt='cheby2_low', output='sos')
x = sc.sosfilt(sos, x)
signal_filt[:, ii] = Gain[ii] * sc.resample(x, n_signal)
## affichage de la réponse du filtre
if plot == True:
w, h = sc.freqz(sos)
plt.semilogx(w/2/np.pi*(fs//n_decimate), 20*np.log10(abs(h)), 'k')
elif n_bande-5 <= ii < n_bande-1:
## de 8 kHz à 16 kHz
sos = oct3_dsgn(fc_oct3[ii], fs, output='sos')
signal_filt[:, ii] = Gain[ii] * sc.sosfilt(sos, signal)
## affichage de la réponse du filtre
if plot == True:
w, h = sc.freqz(sos)
plt.semilogx(w/2/np.pi*fs, 20*np.log10(abs(h)), 'k')
elif n_bande-10 <= ii < n_bande-5:
## de 2,5 kHz à 6,3 kHz
n_decimate = 2
x = sc.decimate(signal, n_decimate)
sos = oct3_dsgn(fc_oct3[ii], fs//n_decimate, output='sos')
x = sc.sosfilt(sos, x)
signal_filt[:, ii] = Gain[ii] * sc.resample(x, n_signal)
## affichage de la réponse du filtre
if plot == True:
w, h = sc.freqz(sos)
plt.semilogx(w/2/np.pi*(fs//n_decimate), 20*np.log10(abs(h)), 'k')
elif n_bande-15 <= ii < n_bande-10:
## de 800 Hz à 2 kHz
n_decimate = 4#round(fs/(2*fmax))-1
x = sc.decimate(signal, n_decimate)
sos = oct3_dsgn(fc_oct3[ii], fs//n_decimate, output='sos')
x = sc.sosfilt(sos, x)
signal_filt[:, ii] = Gain[ii] * sc.resample(x, n_signal)
## affichage de la réponse du filtre
if plot == True:
w, h = sc.freqz(sos)
plt.semilogx(w/2/np.pi*(fs//n_decimate), 20*np.log10(abs(h)), 'k')
elif n_bande-20 <= ii < n_bande-15:
## de 250 Hz à 630 Hz
n_decimate = 8#round(fs/(2*fmax))-1
x = sc.decimate(signal, n_decimate)
sos = oct3_dsgn(fc_oct3[ii], fs//n_decimate, output='sos')
x = sc.sosfilt(sos, x)
signal_filt[:, ii] = Gain[ii] * sc.resample(x, n_signal)
### affichage de la réponse du filtre
if plot == True:
w, h = sc.freqz(sos)
plt.semilogx(w/2/np.pi*(fs//n_decimate), 20*np.log10(abs(h)), 'k')
elif n_bande-25 <= ii < n_bande-20:
## de 80 Hz à 200 Hz
n_decimate = 16#round(fs/(2*fmax))-1
x = sc.decimate(signal, n_decimate)
sos = oct3_dsgn(fc_oct3[ii], fs//n_decimate, output='sos')
x = sc.sosfilt(sos, x)
signal_filt[:, ii] = Gain[ii] * sc.resample(x, n_signal)
## affichage de la réponse du filtre
if plot == True:
w, h = sc.freqz(sos)
plt.semilogx(w/2/np.pi*(fs//n_decimate), 20*np.log10(abs(h)), 'k')
elif n_bande-29 <= ii < n_bande-25:
## de 25 Hz à 63 Hz
n_decimate = 32#round(fs/(2*fmax))-1
x = sc.decimate(signal, n_decimate)
sos = oct3_dsgn(fc_oct3[ii], fs//n_decimate, output='sos')
x = sc.sosfilt(sos, x)
signal_filt[:, ii] = Gain[ii] * sc.resample(x, n_signal)
## affichage de la réponse du filtre
if plot == True:
w, h = sc.freqz(sos)
plt.semilogx(w/2/np.pi*(fs//n_decimate), 20*np.log10(abs(h)), 'k')
if plot == True:
plt.grid(which='both', linestyle='-', color='grey')
# plt.xticks([20, 50, 100, 200, 500, 1000, 2000, 5000, 10000, 20000],
# ["20", "50", "100", "200", "500", "1K",
# "2K", "5K", "10K", "20K"])
plt.xlabel('Fréquence [Hz]'), plt.ylabel('Gain [dB]')
plt.title('Réponse en fréquence des filtres 1/3 d\'octaves')
plt.xlim((10, 22e3)), plt.ylim((-5, 1))
plt.show()
#------- Sommation des signaux filtrés pour recomposer le signal d'origine
S = signal_filt.sum(axis=1)
S = S - np.mean(S)
## tuckey_window = sc.tukey(len(S), alpha=0.01)
## S = tuckey_window * S
G = 10**(gain_general/20)
return G * S
解决方案
推荐阅读
- java - 为什么我没有得到输出？谁能帮我？
- wordpress - 发布帖子时如何将php代码作为内容传递
- python - “元组”对象没有属性“数据”
- python - 向线性程序纸浆添加惩罚参数
- powershell - 如何在powershell脚本中将webrequest与processid相关联
- oracle - 我们可以在 SpringBoot 中使用 Logstash 从 RDBMS 同步数据吗
- swiftui - SwiftUI View 添加 @Environment(\.presentationMode) 崩溃 .sheet() 演示
- scala - Scala Spark 函数，如 group by，describe() 返回不正确的结果
- c++ - 在 C++ 中将类定义拆分为多个模块单元
- javascript - 在 setInterval 循环中使用 setter 函数时，clearInterval 不起作用