Phase shift of the signal after using FIR filter firwin in Python

Question

So after my two last questions I come to my actual problem. Maybe somebody finds the error in my theoretical procedure or I did something wrong in programming.

I am implementing a bandpass filter in Python using scipy.signal (using the firwin function). My original signal consists of two frequencies (w_1=600Hz, w_2=800Hz). There might be a lot more frequencies that's why I need a bandpass filter.

In this case I want to filter the frequency band around 600 Hz, so I took 600 +/- 20Hz as cutoff frequencies. When I implemented the filter and reproduce the signal in the time domain using lfilter the right frequency is filtered.

To get rid of the phase shift I plotted the frequency response by using scipy.signal.freqz with the return h of firwin as numerator and 1 as predefined denumerator. As described in the documentation of freqz I plotted the phase (== angle in the doc) as well and was able to look at the frequency response plot to get the phase shift for the frequency 600 Hz of the filtered signal.

So the phase delay t_p is

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

Unfortunately when I add this phase delay to the time data of my filtered signal, it has not got the same phase as the original 600 Hz signal.

I added the code. It is weird, before eliminating some part of the code to keep the minimum, the filtered signal started at the correct amplitude - now it is even worse.

################################################################################
#
# 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

The delay of your FIR filter is simply 0.5*(n - 1)/fs, where n is the number of filter coefficients (i.e. "taps") and fs is the sample rate. Your implementation of this delay is fine.

The problem is that your array of time values tps is not correct. Take a look at 1.0/(tps[1] - tps[0]); you'll see that it does not equal fs.

Change this:

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

to, for example, this:

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

and your plots of the original and filtered 60 Hz signals will line up beautifully.

---

For another example, see http://wiki.scipy.org/Cookbook/FIRFilter. In the script there, the delay is calculated on line 86. Below this, the delay is used to plot the original signal aligned with the filtered signal.

Note: The cookbook example uses scipy.signal.lfilter to apply the filter. A more efficient approach is to use numpy.convolve.

Answer 2

Seems like you may have had this answered already, but I believe that this is what the filtfilt function is used for. Basically, it does both a forward sweep and a backward sweep through your data, thus reversing the phase shift introduced by the initial filtering. Might be worth looking into.