Low-Pass Filtering of not evenly sampled signal

Question

I have a signal sampled unevenly over 1 million ns. the signal is sampled over 1GHZ clock and the samples are as the following:

0-100 ns - sample every 1 ns.

100-1000 ns - sample every 10 ns.

1000-10,000 ns - every 100 ns.

10,000-100,000 - every 1000 ns.

100,000-1,000,000 - every 10,000 ns.

picture of the sampled (and noisy) signal:

the main goal is to get Z(f) out of the measurments, using the formula:

i'm required first to Filter the signal through a LPF Filter with cutoff frequency of 200 MHz.

for this purpose I decided to try the following: given Time,Voltage measurments:

1. break the data to the evenly sampled segments.

2. for each segment - find sampling frequency.

3. filter each segment seperately and return filtered results.
4. stitch the filtered results to (hopefully) get a full filtered measurments in time.

For this purpose i wrote the following code:

    from __future__ import division
    import numpy as np
    import matplotlib.pyplot as plt
    import numpy.fft as fft
    from scipy import fftpack
    from scipy import interpolate
    from scipy.signal import butter, lfilter, freqz
    import pandas as pd


    def splitToDecades(lTime,lMeasurmentsList,time_units):
        lDecades=[]
        lMeasurmentsPerDecades=[]
        lTempTime=[]
        lTempMeasurments=[]
        res0=lTime[1]-lTime[0]
        for i,m in zip(range(len(lTime)-1),lMeasurmentsList):
            res = lTime[i + 1] - lTime[i]
            if res==res0:
                lTempTime.append(lTime[i])
                lTempMeasurments.append(m)
            else:
                lTempTime.append(lTime[i])
                lTempMeasurments.append(m)
                res0=res
                print res0
                lDecades.append(lTempTime)
                lMeasurmentsPerDecades.append(lTempMeasurments)
                del lTempTime
                del lTempMeasurments
                lTempTime=[]
                lTempMeasurments=[]
        lDecades.append(lTempTime)
        lTempMeasurments.append(lMeasurmentsList[-1])
        lMeasurmentsPerDecades.append(lTempMeasurments)
        del lTempTime
        del lTempMeasurments
        lSamplingFrequencies = []
        for d in lDecades:
            print d
            Ts=d[2]-d[1]
            lSamplingFrequencies.append(1/((Ts)*time_units))
        return lDecades,lMeasurmentsPerDecades,lSamplingFrequencies

    def butter_lowpass(cutoff, fs, order=5):
        nyq = 0.5 * fs
        normal_cutoff = float(cutoff) / nyq
        b, a = butter(order, normal_cutoff, btype='lowpass', analog=False)
        return b, a

    def butter_lowpass_filter(data, cutoff, fs, order=5):
        b, a = butter_lowpass(cutoff, fs, order=order)
        y = lfilter(b, a, data)
        print y
        return y

    def filter_in_parts(lDecades,lMeasurmentsPerDecades,lSamplingFrequencies,cutoff):
        lFilteredData = []
        for lmeas,fs in zip(lMeasurmentsPerDecades,lSamplingFrequencies):
            lFilteredData.append( butter_lowpass_filter(lmeas,cutoff,fs))
        return lFilteredData


    time_units = 1e-9  # ns
    df = pd.read_csv(csv_name)
    Y = df[['Time', 'Voltage']]

    lMeasurmentsList=[]
    lTime=[]
    for t,v in zip(Y['Time'],Y['Voltage']):
        lMeasurmentsList.append(v)
        lTime.append(t)


    lDecades,lMeasurmentsPerDecades,lSamplingFrequencies = splitToDecades(lTime,lMeasurmentsList,time_units)
    lFilteredData =filter_in_parts(lDecades,lMeasurmentsPerDecades,lSamplingFrequencies,200e6)
    lFiltered_final = []
    for p in lFilteredData:
        for d in p:
            lFiltered_final.append(d)

    plt.plot(lFiltered_final)
    plt.show()

I've tried using a butterworth filter and i'm getting the following error:

Traceback (most recent call last):
  File "<module1>", line 77, in <module>
  File "<module1>", line 61, in filter_in_parts
  File "<module1>", line 53, in butter_lowpass_filter
  File "<module1>", line 49, in butter_lowpass
  File "C:\Python27\lib\site-packages\scipy\signal\filter_design.py", line
  2591, in butter
    output=output, ftype='butter', fs=fs)
  File "C:\Python27\lib\site-packages\scipy\signal\filter_design.py", line
  2091, in iirfilter
    raise ValueError("Digital filter critical frequencies "
**ValueError: Digital filter critical frequencies must be 0 < Wn < 1**

can someone please advice what's wrong with my code or alternatively suggest a better approach for filtering such a signal?

Answer 1

I think the problem is that you try to filter signal parts where the filter cutoff frequency is higher than the nyquist rate, e.g. for your last segment, where the sampling period $T_s = 10^4 \cdot 10^{-9}$ seconds corresponds to a sampling frequency of $f_s= 10^-4 \cdot 10^9$ Hz $= 10^5$ Hz. This leads to a nyquist frequency of $\frac{f_s}{2}=5 \cdot 10^4$ Hz, which is way below your desired cutoff frequency of $2\cdot 10^8$ Hz.

My suggestion would be to modify your function butter_lowpass_filter() to first check whether the filtering operation at a given cutoff and sampling frequency would be of any effect, that is if $f_c < \frac{f_s}{2}$. Otherwise, you can just skip the filtering, and the design routine will not be called with faulty parameters (the error appears when your design frequency is not within the appropriate range).

Also, I am not sure whether you can just filter individual segments like this, since the effects from smoothing a segment might also appear in successive segments. Using your code, this is not possible. I would suggest doing interpolation of the segments with lower sampling rates first, for example by transforming them into frequency domain, appending zeros there, and then transforming them back into the time domain. That way, you are not adding any information, but getting more samples which would enable you to filter the signal at uniform sampling rate.

The number of additional zeros is different for each segment. However, you may be able to use the scipy.signal.resample()-routine which i believe does all of that for you.

For example, the second segment has a sampling rate of $f_{s,2} =\frac{f_{s,1}}{10}$, so it has to be upsampled to the same sampling rate of the first segment. For the second segment this means there will be ten times as much samples. For the rest of the segments, this will vary according to their sampling rates, but you can write a function to determine the number of samples:

$$N_\text{os} = N_\text{original} \cdot \frac{f_{s,1}}{f_{s, l}},$$ where $l$ is the segment number, and $f_{s,1}$ corresponds to the sampling rate of the first segment (the one with the highest sampling frequency) $N_{os}$ denotes the number of samples that should be the result from oversampling, and $N_\text{original}$ the original number of samples in the segment $l$.