在Python中使用FIR滤波器firwin后的信号相移

Question

所以在我最后两个问题后，我来到我的实际问题。也许有人在我的理论程序中发现错误，或者我在编程时做了错误。

我正在使用scipy.signal（使用firwin函数）在Python中实现带通滤波器。我的原始信号由两个频率组成（w_1 = 600Hz，w_2 = 800Hz）。可能有更多的频率，这就是为什么我需要一个带通滤波器。

在这种情况下，我想过滤600赫兹左右的频带，所以我把600 +/- 20Hz作为截止频率。当我使用lfilter实现滤波器并在时域中重现信号时，正确的频率被滤波。

为了摆脱相移，我使用scipy.signal.freqz绘制频率响应，firwin的返回值h为分子，1为预定义的分子。 正如freqz文档中所描述的那样，我也绘制了相位（文档中的==角度），并能够查看频率响应图，以获得滤波信号频率600 Hz的相移。

所以相位延迟T_P是

T_P = - （Tetha（W））/（W）

不幸的是，当我这个相位延迟添加到我的滤波后的信号的时间数据，它没有与原始的600Hz信号具有相同的相位。

我添加了代码。奇怪的是，在消除代码的一部分以保持最小值之前，滤波后的信号以正确的幅度开始 - 现在情况更糟。

################################################################################ 
# 
# Filtering test 
# 
################################################################################ 
# 
from math import * 
import numpy as np 
from scipy import signal 
from scipy.signal import firwin, lfilter, lti 
from scipy.signal import freqz 

import matplotlib.pyplot as plt 
import matplotlib.colors as colors 

################################################################################ 
# Nb of frequencies in the original signal 
nfrq = 2 
F = [60,80] 

################################################################################ 

# Sampling: 
nitper = 16 
nper = 50. 

fmin = np.min(F) 
fmax = np.max(F) 

T0 = 1./fmin 
dt = 1./fmax/nitper 

#sampling frequency 
fs = 1./dt 
nyq_rate= fs/2 

nitpermin = nitper*fmax/fmin 
Nit = int(nper*nitpermin+1) 
tps = np.linspace(0.,nper*T0,Nit) 
dtf = fs/Nit 

################################################################################ 

# Build analytic signal 
# s = completeSignal(F,Nit,tps) 
scomplete = np.zeros((Nit)) 
omg1 = 2.*pi*F[0] 
omg2 = 2.*pi*F[1] 
scomplete=scomplete+np.sin(omg1*tps)+np.sin(omg2*tps) 

#ssingle = singleSignals(nfrq,F,Nit,tps) 
ssingle=np.zeros((nfrq,Nit)) 
ssingle[0,:]=ssingle[0,:]+np.sin(omg1*tps) 
ssingle[1,:]=ssingle[0,:]+np.sin(omg2*tps) 
################################################################################ 

## Construction of the desired bandpass filter 
lowcut = (60-2) # desired cutoff frequencies 
highcut = (60+2) 
ntaps = 451  # the higher and closer the signal frequencies, the more taps for the filter are required 

taps_hamming = firwin(ntaps,[lowcut/nyq_rate, highcut/nyq_rate], pass_zero=False) 

# Use lfilter to get the filtered signal 
filtered_signal = lfilter(taps_hamming, 1, scomplete) 

# The phase delay of the filtered signal 
delay = ((ntaps-1)/2)/fs 

plt.figure(1, figsize=(12, 9)) 
# Plot the signals 
plt.plot(tps, scomplete,label="Original signal with %s freq" % nfrq) 
plt.plot(tps-delay, filtered_signal,label="Filtered signal %s freq " % F[0]) 
plt.plot(tps, ssingle[0,:],label="original signal %s Hz" % F[0]) 
plt.grid(True) 
plt.legend() 
plt.xlim(0,1) 
plt.xlabel('Time (s)') 
plt.ylabel('Amplitude') 

# Plot the frequency responses of the filter. 
plt.figure(2, figsize=(12, 9)) 
plt.clf() 
# First plot the desired ideal response as a green(ish) rectangle. 
rect = plt.Rectangle((lowcut, 0), highcut - lowcut, 5.0,facecolor="#60ff60", alpha=0.2,label="ideal filter") 
plt.gca().add_patch(rect) 
# actual filter 
w, h = freqz(taps_hamming, 1, worN=1000) 
plt.plot((fs * 0.5/np.pi) * w, abs(h), label="designed rectangular window filter") 

plt.xlim(0,2*F[1]) 
plt.ylim(0, 1) 
plt.grid(True) 
plt.legend() 
plt.xlabel('Frequency (Hz)') 
plt.ylabel('Gain') 
plt.title('Frequency response of FIR filter, %d taps' % ntaps) 
plt.show()' 

Answer 1

您的FIR滤波器的延迟仅仅是0.5*(n - 1)/fs，其中n是（即“轻敲”）的滤波器系数的数目和fs是采样率。您的延迟执行情况良好。

问题是您的时间值数组tps不正确。看看 在1.0/(tps[1] - tps[0]);你会发现它并不等于fs。

更改此：

tps = np.linspace(0.,nper*T0,Nit) 

，例如，这样的：

T = Nit/fs 
tps = np.linspace(0., T, Nit, endpoint=False) 

和你原来的情节和滤波60 Hz的信号将精美排队。

---

又如，参见http://wiki.scipy.org/Cookbook/FIRFilter。 在那里的脚本中，延迟是在第86行计算的。在此之下，延迟用于绘制与滤波信号对齐的原始信号。

注意：食谱示例使用scipy.signal.lfilter来应用过滤器。更有效的方法是使用numpy.convolve。

Answer 2

似乎你可能已经得到了这个答案，但我相信这是filtfilt函数用于什么。基本上，它通过数据执行前向扫描和后向扫描，从而逆转由初始滤波引起的相移。可能值得研究。