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

2 投票
2 回答
4896 浏览
提问于 2025-04-18 00:08

在我问完之前的两个问题后,我现在遇到了一个实际的问题。也许有人能发现我理论过程中的错误,或者我在编程时做错了什么。

我正在用Python实现一个带通滤波器,使用的是scipy.signal库中的firwin函数。我的原始信号包含两个频率(w_1=600Hz,w_2=800Hz)。可能还有很多其他频率,这就是我需要带通滤波器的原因。

在这种情况下,我想过滤掉600 Hz附近的频率,所以我选择了600 +/- 20Hz作为截止频率。当我实现滤波器并使用lfilter在时间域中重现信号时,正确的频率被过滤掉了。

为了消除相位偏移,我使用scipy.signal.freqz绘制了频率响应,使用firwin的返回值h作为分子,1作为预定义的分母。根据freqz的文档,我也绘制了相位(文档中的角度),并能够查看频率响应图,以获取过滤信号在600 Hz频率下的相位偏移。

因此,相位延迟t_p是

t_p=-(Tetha(w))/(w)

不幸的是,当我将这个相位延迟添加到过滤信号的时间数据中时,它的相位与原始的600 Hz信号并不相同。

我添加了代码。很奇怪,在删除一些代码以保持最小化之前,过滤信号的幅度是正确的——现在情况甚至更糟。

################################################################################
#
# 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()'

2 个回答

2

看起来你可能已经得到了这个问题的答案,但我觉得这就是 filtfilt 函数的用途。简单来说,它会对你的数据进行一次正向处理和一次反向处理,这样就能消除最初过滤时产生的相位偏移。值得你去了解一下。

2

你的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赫兹信号就能完美对齐了。


想要更多例子,可以看看这个链接:http://wiki.scipy.org/Cookbook/FIRFilter。在那里的脚本中,延迟是在第86行计算的。接下来,这个延迟被用来将原始信号和过滤后的信号对齐进行绘图。

注意:这个食谱示例使用scipy.signal.lfilter来应用滤波器。更高效的方法是使用numpy.convolve

撰写回答