【问题标题】:Phase shift of the signal after using FIR filter firwin in Python在 Python 中使用 FIR 滤波器 firwin 后信号的相移
【发布时间】:2014-05-06 10:11:31
【问题描述】:

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

我正在使用scipy.signal(使用 firwin 函数)在 Python 中实现带通滤波器。我的原始信号由两个频率组成(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()'

【问题讨论】:

    标签: python scipy filtering


    【解决方案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

    【讨论】:

    • 谢谢!!!实际上我在里面已经有这个延迟,但它不起作用。我认为这些是损坏的样本,并且这个 t_p 延迟是额外的。我在 Sundararajan,信号处理中找到了这个。
    • “它不起作用”是什么意思?如果您可以编辑您的问题以添加您尝试过的代码、您得到的结果以及解释为什么它不是您所期望的,这将有所帮助。
    • 我更新了我的答案。感谢您将代码添加到问题中。这使得发现问题变得更加容易。
    • 你拯救了我的周末!!!我认为错误在一个完全不同的地方。尽管如此,我想知道为什么相移仅取决于标签大小,而不取决于零极点图中零点的位置(因此是频率响应)。我也试试卷积函数。
    【解决方案2】:

    似乎您可能已经回答了这个问题,但我相信这就是 filtfilt 函数的用途。基本上,它对您的数据进行前向扫描和后向扫描,从而逆转初始滤波引入的相移。可能值得研究。

    【讨论】:

      猜你喜欢
      • 2014-04-24
      • 1970-01-01
      • 2014-04-26
      • 1970-01-01
      • 1970-01-01
      • 1970-01-01
      • 1970-01-01
      • 1970-01-01
      • 1970-01-01
      相关资源
      最近更新 更多