Signal Processing Python

[Malva Ferster]

Thesignal processing toolbox currently contains some filteringfunctions, a limited set of filter design tools, and a few B-splineinterpolation algorithms for 1- and 2-D data. While theB-spline algorithms could technically be placed under theinterpolation category, they are included here because they only workwith equally-spaced data and make heavy use of filter-theory andtransfer-function formalism to provide a fast B-spline transform. Tounderstand this section, you will need to understand that a signal inSciPy is an array of real or complex numbers.

A B-spline is an approximation of a continuous function over a finite-domain in terms of B-spline coefficients and knot points. If the knot-points are equally spaced with spacing \(\Delta x\), then the B-splineapproximation to a 1-D function is the finite-basis expansion.

In these expressions, \(\beta^o\left(\cdot\right)\) is the space-limitedB-spline basis function of order \(o\). The requirement of equally-spacedknot-points and equally-spaced data points, allows the development of fast(inverse-filtering) algorithms for determining the coefficients, \(c_j\),from sample-values, \(y_n\). Unlike the general spline interpolationalgorithms, these algorithms can quickly find the spline coefficients for largeimages.

The advantage of representing a set of samples via B-spline basisfunctions is that continuous-domain operators (derivatives, re-sampling, integral, etc.), which assume that the data samples are drawnfrom an underlying continuous function, can be computed with relativeease from the spline coefficients. For example, the second derivativeof a spline is

The savvy reader will have already noticed that the data samples are relatedto the knot coefficients via a convolution operator, so that simpleconvolution with the sampled B-spline function recovers the original data fromthe spline coefficients. The output of convolutions can change depending onhow the boundaries are handled (this becomes increasingly more important as thenumber of dimensions in the dataset increases). The algorithms relating toB-splines in the signal-processing subpackage assume mirror-symmetricboundary conditions. Thus, spline coefficients are computed based on thatassumption, and data-samples can be recovered exactly from the splinecoefficients by assuming them to be mirror-symmetric also.

Filtering is a generic name for any system that modifies an inputsignal in some way. In SciPy, a signal can be thought of as a NumPyarray. There are different kinds of filters for different kinds ofoperations. There are two broad kinds of filtering operations: linearand non-linear. Linear filters can always be reduced to multiplicationof the flattened NumPy array by an appropriate matrix resulting inanother flattened NumPy array. Of course, this is not usually the bestway to compute the filter, as the matrices and vectors involved may behuge. For example, filtering a \(512 \times 512\) image with thismethod would require multiplication of a \(512^2 \times 512^2\)matrix with a \(512^2\) vector. Just trying to store the\(512^2 \times 512^2\) matrix using a standard NumPy array wouldrequire \(68,719,476,736\) elements. At 4 bytes per element thiswould require \(256\textrmGB\) of memory. In most applications,most of the elements of this matrix are zero and a different methodfor computing the output of the filter is employed.

Many linear filters also have the property of shift-invariance. Thismeans that the filtering operation is the same at different locationsin the signal and it implies that the filtering matrix can beconstructed from knowledge of one row (or column) of the matrix alone.In this case, the matrix multiplication can be accomplished usingFourier transforms.

This same function convolve can actually take N-Darrays as inputs and will return the N-D convolution of thetwo arrays, as is shown in the code example below. The same input flags areavailable for that case as well.

Calculating the convolution in the time domain as above is mainly used forfiltering when one of the signals is much smaller than the other ( \(K\ggM\) ), otherwise linear filtering is more efficiently calculated in thefrequency domain provided by the function fftconvolve. By default,convolve estimates the fastest method using choose_conv_method.

The difference-equation filter is called using the command lfilter inSciPy. This command takes as inputs the vector \(b,\) the vector,\(a,\) a signal \(x\) and returns the vector \(y\) (the samelength as \(x\) ) computed using the equation given above. If \(x\) isN-D, then the filter is computed along the axis provided.If desired, initial conditions providing the values of\(z_0\left[-1\right]\) to \(z_K-1\left[-1\right]\) can be providedor else it will be assumed that they are all zero. If initial conditions areprovided, then the final conditions on the intermediate variables are alsoreturned. These could be used, for example, to restart the calculation in thesame state.

Sometimes, it is more convenient to express the initial conditions in terms ofthe signals \(x\left[n\right]\) and \(y\left[n\right].\) In otherwords, perhaps you have the values of \(x\left[-M\right]\) to\(x\left[-1\right]\) and the values of \(y\left[-N\right]\) to\(y\left[-1\right]\) and would like to determine what values of\(z_m\left[-1\right]\) should be delivered as initial conditions to thedifference-equation filter. It is not difficult to show that, for \(0\leqm

Time-discrete filters can be classified into finite response (FIR) filters andinfinite response (IIR) filters. FIR filters can provide a linear phaseresponse, whereas IIR filters cannot. SciPy provides functionsfor designing both types of filters.

Note that firwin uses, per default, a normalized frequency defined suchthat the value \(1\) corresponds to the Nyquist frequency, whereas thefunction freqz is defined such that the value \(\pi\) correspondsto the Nyquist frequency.

SciPy provides two functions to directly design IIR iirdesign andiirfilter, where the filter type (e.g., elliptic) is passed as anargument and several more filter design functions for specific filter types,e.g., ellip.

The example below designs an elliptic low-pass filter with defined pass-bandand stop-band ripple, respectively. Note the much lower filter order (order 4)compared with the FIR filters from the examples above in order to reach the samestop-band attenuation of \(\approx 60\) dB.

The sos format is a single 2-D array of shape (n_sections, 6),representing a sequence of second-order transfer functions which, whencascaded in series, realize a higher-order filter with minimal numericalerror. Each row corresponds to a second-order tf representation, withthe first three columns providing the numerator coefficients and the lastthree providing the denominator coefficients:

The coefficients are typically normalized, such that \(a_0\) is always 1.The section order is usually not important with floating-point computation;the filter output will be the same, regardless of the order.

The IIR filter design functions first generate a prototype analog low-pass filterwith a normalized cutoff frequency of 1 rad/sec. This is then transformed intoother frequencies and band types using the following substitutions:

A median filter is commonly applied when noise is markedly non-Gaussian orwhen it is desired to preserve edges. The median filter works by sorting allof the array pixel values in a rectangular region surrounding the point ofinterest. The sample median of this list of neighborhood pixel values is usedas the value for the output array. The sample median is the middle-array valuein a sorted list of neighborhood values. If there are an even number ofelements in the neighborhood, then the average of the middle two values isused as the median. A general purpose median filter that works onN-D arrays is medfilt. A specialized version that worksonly for 2-D arrays is available as medfilt2d.

A median filter is a specific example of a more general class of filterscalled order filters. To compute the output at a particular pixel, all orderfilters use the array values in a region surrounding that pixel. These arrayvalues are sorted and then one of them is selected as the output value. Forthe median filter, the sample median of the list of array values is used asthe output. A general-order filter allows the user to select which of thesorted values will be used as the output. So, for example, one could choose topick the maximum in the list or the minimum. The order filter takes anadditional argument besides the input array and the region mask that specifieswhich of the elements in the sorted list of neighbor array values should beused as the output. The command to perform an order filter isorder_filter.

The Wiener filter is a simple deblurring filter for denoising images. This isnot the Wiener filter commonly described in image-reconstruction problems but,instead, it is a simple, local-mean filter. Let \(x\) be the input signal,then the output is

where \(m_x\) is the local estimate of the mean and\(\sigma_x^2\) is the local estimate of the variance. The window forthese estimates is an optional input parameter (default is \(3\times3\) ).The parameter \(\sigma^2\) is a threshold noise parameter. If\(\sigma\) is not given, then it is estimated as the average of the localvariances.

Spectral analysis refers to investigating the Fourier transform [1] of asignal. Depending on the context, various names, like spectrum, spectral density orperiodogram exist for the various spectral representations of the Fourier transform.[2] This section illustrates the most common representations by the exampleof a continuous-time sine wave signal of fixed duration. Then the use of the discreteFourier transform [3] on a sampled version of that sine wave is discussed.

Note that the concept of Fourier series is closely related but differs in a crucialpoint: Fourier series have a spectrum made up of discrete-frequency harmonics, while inthis section the spectra are continuous in frequency.

3a8082e126