2013. január 17., csütörtök

Implementing a Hilbert (90 degree shift) filter in Python

Why?

A digital 90° phase shift filter is an important building block of the so-called Software Defined Radios (SDRs).

And SDR is a radio that has (relatively) minimal hardware, and most of the features are implemented in software. There are no fancy buttons and displays, but there is a UI application that controls the box. The box is a direct conversion receiver. It just converts part of the radio frequency spectrum to the audio spectrum (the reality is a bit more complicated, but you got the idea).

An SDR usually provides two signals that are almost the same, except one signal's frequency components are shifted 90 degrees in one direction.

Given these two audio signals, the software can do anything a conventional receiver can achieve with bulky components or expensive integrated circuits. Plus, the software can be changed easily. Actually this is the greatest thing in SDRs.

The software can demodulate AM signals by LW-MW-SW broadcast stations, or FM usually encountered in the URH spectrum. It can make amateur SSB and CW signals readable, or can decode digital transmissions like DAB. Anything you can imagine.



SSB reception

I'm mainly interested in amateur (HAM) radio, and so I'd like to receive SSB transmissions. To do that, I have to 90° phase-shift one of the incoming audio signals from the SDR box, and sum with the other. The result is an SSB receiver.

FIR filters

I use the numpy.filter Python package to process audio signals. The firwin and firwin2 function are very useful for designing all sorts of FIR filters, but I could not find a built-in function that can readily be used to shift all frequencies by 90 degrees.

A filter is represented by its coefficients. The simplest FIR filter one can imagine is the averaging filter. You take N samples, add them, and divide the result by N. If you encounter a new sample, you just drop one sample from one end of the line, and put the new sample to the other end. You add them and divide by N again. Simple. The general FIR filter is a weighted averaging filter. The coefficients are the weights.

The coefficients are also the filter's impulse response. That is, if you filter a signal that has all zero samples except one, you get all zeroes and the coefficients somewhere in the middle. Example:

Signal:
0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 ...

Coefficients:
0.3 0.4 0.3

Response:

0 0 0 0 0 0 0 0 0 0.3 0.4 0.3 0 ...

As you can see, the response contains the coefficients.

This is important because the impulse response can be calculated from the desired frequency response by the inverse discrete Fourier Transform.

90° phase shift

 In our case, this leads us to the code below:

import numpy.fft
import matplotlib.pyplot as plt

coeffs = numpy.fft.fftshift(
    numpy.fft.ifft([0]+[1]*100+[0]*100)
)

plt.plot(numpy.imag(coeffs))
plt.plot(numpy.real(coeffs))
plt.show()

The variable coeffs holds the filter coefficients. It can be used with the numpy.signal.lfilter function to process data.

The code works by calculating the inverse discrete Fourier Transform of a strange frequency response. It reads like this: "pass every negative frequencies, and supress all of the positive frequencies". This series has a complex iDFT. We are interested in the imaginary part of this inverse. The fftshift function rearranges the result. The plt.* function calls shows us the result.

Note that the filter coefficients are complex. An audio signal from and SDR receiver can be processed by taking the left and rigth channel as the real and imaginary parts of a complex signal, and filtering that - also complex - signal. The filtered signal will also be complex, and the imaginary part will be shifted by 90° (at every frequency).

If the signal to be filtered is real-valued, then the resulting filtered signal will still be a complex one: a real, and a 90°-shifted imaginary.