«

»

Jul 09

Nyquist Signal Expansion with Python

Filtering and Scaling

The real magic is the finite impulse response filter, or FIR filter.  We are going to cheat and let the computer figure out what the filter coefficients need to be.  To do this, we use the remez() function in the scipy.signals package.  This function is similar to the firpm() function in MATLAB.  We give the function our pass band and stop band.  We want to pass all frequencies less than or equal to fmax.  Since we expanded the signal, we need to compute what fmax is for the higher sample rate.  To do this we divide by L.  This gives us:

h = L*signal.remez(20*L+1, [0, 2/float(fs_in*L), 4/float(fs_in*L), 0.5], [1, 0],[1, 1])

The magnitude of the filter is:

Magnitude response of generated filter.

Figure 3.  Magnitude response of generated filter.

If you noticed I snuck a scaling factor into the filter coefficients.  I scaled them by L.

Now we need to run the expanded signal through the filter:

x_ef=signal.lfilter(h,1,x_e)

Which gives us:

Expanded data run through the filter with original signal in background.

Figure 4.  Expanded data run through the filter with original signal in background.

Now, this looks like it could be what we are looking for, but why is it shifted you ask?  That has to do with the delay in the filter, which is directly related to the number of taps in the filter.

Calculating the Filter Delay

The filter delay is equal to the number of taps, minus 1, in the filter divided by 2.  For this example it is n_s=\frac{\textrm{len}(h)-1}{2}.  That’s the delay in samples, but how does that relate back to the time shift? 

delay = Ts_in/float(L)*(len(h)-1)/2.0

.  Now that we’ve figured this out, we will shift the input signal to match the delay coming out of the filter:

Filtered, expanded data with the initial signal delayed to match.

Figure 5.  Filtered, expanded data with the initial signal delayed to match.

The calculated delay for this example is 1.66.  Looking at Figure 5, we see that the points match up after that time.

Non-integer Resampling

So you say, this is great, but I need my signal resampled by a non-integer value.  That’s great.  We do the exact same thing as above and then slice the resulting data to take every Mth sample.  Easy as

x_sliced = x_ef[0::M]

Next we put it all together.

Permanent link to this article: http://blog.curioussystem.com/2013/07/nyquist-signal-expansion-with-python/

Leave a Reply

Your email address will not be published. Required fields are marked *

You may use these HTML tags and attributes: <a href="" title=""> <abbr title=""> <acronym title=""> <b> <blockquote cite=""> <cite> <code> <del datetime=""> <em> <i> <q cite=""> <s> <strike> <strong>