parameters for low pass fir filter using scipy

Question

I am trying to write a simple low pass filter using scipy, but I need help defining the parameters.

I have 3.5 million records in the time series data that needs to be filtered, and the data is sampled at 1000 hz.

I am using signal.firwin and signal.lfilter from the scipy library.

The parameters I am choosing in the code below do not filter my data at all. Instead, the code below simply produces something that graphically looks like the same exact data except for a time phase distortion that shifts the graph to the right by slightly less than 1000 data points (1 second).

In another software program, running a low pass fir filter through graphical user interface commands produces output that has similar means for each 10 second (10,000 data point) segment, but that has drastically lower standard deviations so that we essentially lose the noise in this particular data file and replace it with something that retains the mean value while showing longer term trends that are not polluted by higher frequency noise. The other software's parameters dialog box contains a check box that allows you to select the number of coefficients so that it "optimizes based on sample size and sampling frequency." (Mine are 3.5 million samples collected at 1000 hz, but I would like a function that uses these inputs as variables.)

*Can anyone show me how to adjust the code below so that it removes all frequencies above 0.05 hz?* I would like to see smooth waves in the graph rather than just the time distortion of the same identical graph that I am getting from the code below now.

class FilterTheZ0():
    def __init__(self,ZSmoothedPylab):
        #------------------------------------------------------
        # Set the order and cutoff of the filter
        #------------------------------------------------------
        self.n = 1000
        self.ZSmoothedPylab=ZSmoothedPylab
        self.l = len(ZSmoothedPylab)
        self.x = arange(0,self.l)
        self.cutoffFreq = 0.05

        #------------------------------------------------------
        # Run the filter
        #------------------------------------------------------
        self.RunLowPassFIR_Filter(self.ZSmoothedPylab, self.n, self.l
                                       , self.x, self.cutoffFreq)

    def RunLowPassFIR_Filter(self,data, order, l, x, cutoffFreq):
        #------------------------------------------------------
        # Set a to be the denominator coefficient vector
        #------------------------------------------------------
        a = 1
        #----------------------------------------------------
        # Create the low pass FIR filter
        #----------------------------------------------------
        b = signal.firwin(self.n, cutoff = self.cutoffFreq, window = "hamming")

        #---------------------------------------------------
        # Run the same data set through each of the various
        # filters that were created above.
        #---------------------------------------------------
        response = signal.lfilter(b,a,data)
        responsePylab=p.array(response)

        #--------------------------------------------------
        # Plot the input and the various outputs that are
        # produced by running each of the various filters
        # on the same inputs.
        #--------------------------------------------------

        plot(x[10000:20000],data[10000:20000])
        plot(x[10000:20000],responsePylab[10000:20000])
        show()

        return

Answer 1

Cutoff is normalized to the Nyquist frequency, which is half the sampling rate. So with FS = 1000 and FC = 0.05, you want cutoff = 0.05/500 = 1e-4.

from scipy import signal

FS = 1000.0                                          # sampling rate
FC = 0.05/(0.5*FS)                                   # cutoff frequency at 0.05 Hz
N = 1001                                             # number of filter taps
a = 1                                                # filter denominator
b = signal.firwin(N, cutoff=FC, window='hamming')    # filter numerator

M = FS*60                                            # number of samples (60 seconds)
n = arange(M)                                        # time index
x1 = cos(2*pi*n*0.025/FS)                            # signal at 0.025 Hz
x = x1 + 2*rand(M)                                   # signal + noise
y = signal.lfilter(b, a, x)                          # filtered output

plot(n/FS, x); plot(n/FS, y, 'r')                    # output in red
grid()

The filter output is delayed half a second (the filter is centered on tap 500). Note that the DC offset added by the noise is preserved by the low-pass filter. Also, 0.025 Hz is well within the pass range, so the output swing from peak to peak is approximately 2.