Andrew Que Sites list Photos
Projects Contact
Main

November 06, 2017

Impulse Noise Filtering - Implementation of Median Filter

When implementing the median filter there are two options: off-line or real-time. If not done off-line then the filter loops though each data point and calculates the median over a window of the data. Here is an example with a buffer size of 5 and the median window size of 3:

55

42

1000

86

61

55

42

1000

86

61

55

42

1000

86

61

Taking the median for the number highlighted in blue for each row, the filtered data would then be:

55

86

86

Note how the spike of 1000 is filtered out. Note also the resulting buffer is shorter. The filtered data size is always the initial data size less the window size.

If running the filter in real-time, there needs to be a ring-buffer that holds the window size samples. The filtered result is the median of this buffer. New samples replace old samples and the filter continues. Starting this filter requires either loading some starting values or waiting until enough samples have been taken to fill the buffer.

If the following is the initial ring buffer at start:

21

30

71

62

72

Then the first output is 62. Assume the current into into the buffer is the second sample. When a new sample arrives, it will replace the second sample. For instance if the new sample is 41, the ring buffer becomes:

21

41

71

62

72

The median of this is again 62 which is the output of the filter. In this way the filter can be used in real-time.

Getting the filter started is up to the implementer. Most of the time simply using an initial buffer full of zeros works just fine. One can also delay producing filtered output until enough samples have been acquired. The entire buffer can be loaded entirely with the first sample taken. Or the window used for calculations can gradually grow until it reaches the desired size. Which option used depends on the application.

Basic implementation of the median calculation itself is simple. Sort the data and choose the middle data point. However sorting data is a time consuming task and fast sorting algorithms don’t have fixed run times. In addition if one is using a ring-buffer to hold samples the sorting cannot be done in place.

The easiest filter to implement is a selection sort. It is simple and always takes the same number of steps for any window size. Selection sort is also the slowest sort with a fixed complexity of O( n2 ). One bit of acceleration can be gained by only sorting half the data. Since the median is looking for the middle data point in a sorted set, one needs only to find the middle data point—everything else is ignored.

My co-worker Flux suggested I look into using an insertion sort for small filter window sizes. Worst case an insertion sort has the same speed, O( n2 ) as a selection sort. However the speed improves the more sorted the data initially is. If already sorted the speed of the algorithm is O( n ). Since we are running as a windowed filter we are only ever changing one piece of data at a time. Thus the data remains mostly sorted. As a result the insertion sort runs closer to O( n ) complexity and worst case O( 2n )—once through the array, and once again needing to move the last element to the first element. That allows this filter to still be useful in real-time applications, especially for small window sizes.

Here is a simple C implementation of a median filter useful for real-time applications.

//-----------------------------------------------------------------------------
// Uses:
//   Return the median of the data buffer.
// Input:
//   data - Data from which to select median.
//   indexes - Buffer to hold sort indexes (must contain 'windowSize' elements).
//   windowSize - Number of data points.
// Output:
//   Median of data points.
// Notes:
//   Designed for small data sizes.
//   The 'indexes' must be pre-initialized at least once.  After
//   initialization it can be used for all subsequent calls of the same
//   window size.
//       for ( int index = 0; index < windowSize; index += 1 )
//         indexes[ index ] = index;
// Author:
//   Andrew Que <https://www.DrQue.net/>
//-----------------------------------------------------------------------------
static inline int median
(
  int const * data,
  int * indexes,
  int windowSize
)
{
  // Insert sort the data using swap buffer.
  for ( int indexA = 1; indexA < windowSize; indexA += 1 )
  {
    int value = data[ indexes[ indexA ] ];
    int swapIndex = indexes[ indexA ];

    int indexB = indexA - 1;
    while ( ( indexB >= 0 )
         && ( data[ indexes[ indexB ] ] > value ) )
    {
      indexes[ indexB + 1 ] = indexes[ indexB ];
      indexB -= 1;
    }

    indexes[ indexB + 1 ] = swapIndex;
  }

  // Return median.
  return data[ indexes[ windowSize >> 1 ] ];
}

The code is designed for a ring-buffer and requires a second buffer to hold indexes. One should replace all int with the necessary word size for the application but otherwise the code should be drop-in. The algorithm is fairly efficient up to a few tens of window size. After that some more efficient sorting techniques should be used.

November 05, 2017

Impulse Noise Filtering - Median Filter Demo

 
 

In this demonstration there are three sliders. From left to right they are: signal noise, spike density, and filter window length. The signal noise is just basic white noise added to the signal. The spike density is the percentage of impulse noise added to the signal. And the filter window is how many samples are used by the median filter.

Shown in red is the raw signal, complete with spikes. In blue is the median filtered signal. The goal is to have the blue signal from a smooth trace with no spikes, and the default settings should produce this. Increasing the spike density will require a larger filter window.

Just as with an average filter, the median filter will induce phase shift as well as attenuate the signal. The larger the filter window the larger the phase shift and the more attenuation. Unique to the median filter is the creation of plateaus around signal peaks. This is because of the phase shift and the changing direction of the signal lags during peak values.

November 04, 2017

Impulse Noise Filtering - White Noise vs Impulse Noise

There are two main types of noise one encounters in the embedded software world. They are white noise and impulse noise. Let’s start with a perfect signal.

Here is a plain sine wave with no distortion. For all examples, this is the actual signal. We will start with white noise. This is a mainly uniform random addition to the signal of varying amplitude.

Here we see the same signal but with the addition of white noise. Some amount of white noise is common in analog signals. Often a hardware low-pass filter is added to help remove it, but white noise can be the result of the analog to digital converter. Filtering such noise is easy with a software low-pass filter and I have written in the past about using a rolling average to do such filtering.

This is a simple windowed average that has reconstructed the original sine wave. Note that because of the window size there is a lag in the beginning. Also the entire signal has been shifted (phase shift) to the right and the amplitude has been attenuated. These are artifacts of using this kind of low-pass filter. Typically these artifacts do not present too much of a problem but it really depends on the signal being filtered.

There are many sources of white noise and good hardware practices often reduce much of this noise before it ever gets to an analog to digital converter. Nonetheless this is the most common noise I encounter in embedded projects.

Another kind of noise electrical impulse noise. These are random rapid changes in the signal often called spikes. Unlike white noise the magnitude is more dramatic.

This is an example of impulse noise. Filtering noise impulse noise is more difficult because the simple average low-pass filter isn’t as effective at removing this kind of noise.

Shown is a windowed-average low-pass filter being applied to the signal with impulse noise. Averages are thrown off by large deviations and thus the filtered signal is fairly messy.

One solution is to instead use a windowed median filter. It is a low-pass filter but using the median rather than the average. The median is simply the middle value of a set of sorted numbers and it is trivial to calculate. For example in the following set of number: 34, 96, 63, 74, 81. First sort the list: 34, 63, 74, 81, 96. Then pick the middle value which is 74.

Why this is useful in an impulse noise is because a large spike in values will always end up at the beginning or end of the sorted list, and because the median always chooses from the middle these outlier are always ignored. If we were to use a standard first-order low-pass filter (i.e. the average) this would not be the case. For example, say we have the following numbers: 34, 81, 63, 1000, 82. Clearly the large spike is the value 1000, which when sorted ends up as the highest value. The median of this set is 81 where as the average is 252.

This is a median filter in orange with the average filter in green. The median filter does a far better job of reconstructing the original signal. Median filters also introduce phase shift and signal attenuation. In addition they produce plateaus around local extrema.

This signal has a large filter window and has strong phase shift to the right, some signal attenuation, and the peaks are all plateaued. Again, these artifacts need to be understood before using a median filter. The smaller the filter window size, the less the effect. However, if the filter window is too small noise spikes can get through.

Here the window size was too small and while many of the spikes have been removed, some have not.

In subsequent articles I will cover how to implement a median filter and give a demonstration of how it works.

   I ordered a new jacket a couple days ago hoping to find something that would work better for biking in the 30 degree temperature range.  At 38°F I can ride using my heaviest cycling shirt, but I know it won't be good enough much colder than that.  With winter on its way it will be much colder than that.  I started in the morning with the jacket and a thin long sleeve shirt.  Fairly early into the ride I could tell I was going to be too hot.  About 2 miles in I opened the front of the jacket to cool off.  For the ride home I wore just the jacket and no undershirt.  That didn't seem to help either. 
   The jacket is a nice wind break, and is quite warm.  It will be able to handle cold days just fine, but the upper 30s are too warm.  I will have to search for another solution.
   Pictured is sunrise as seen from the west shore of Lake Monona.
   Second year of turning the garage into a Mad Science Laboratory and this Halloween we had a great turnout.  As usual I didn't get everything in the set like I wanted, but what we did have went over well.  Everyone who came though loved it.
   Pictured is Dr. Xiphos trying to sedate the monster.  If he breaks out of the chains again there is no telling who he's gonna eat this time!
   Took a chilly bike ride into work and an even chillier ride home.  The temperature both times was 38°F (3.3°C) but I had a very strong headwind on the ride home.  The wind was high enough (14 MPH with guests of 22 MPH) my average speed dropped from around 14 MPH to around 12.5 MPH (around 10% drop).  On the way back it even started to sleet and was the first freezing precipitation I've encountered this season.  Still, I managed to do with the ride home wearing just a long sleeve shirt, fitted pants, and gloves.  While it was rather cold getting started, once my heart rate was up I was as comfortable as one can be in winds like that.  Yet I was pretty pleased to be riding.
   The ride home would complete updated monthly goal, and on arriving home I finished 605.3 miles for the month of October.  I do not plan to ride tomorrow so this will be the finial for the month.  That easily beasts my previous record set last month of 497.3 miles (which I initially thought was 500 miles).  I rode a total of 46 hours and 49 minutes, burned around and 33,500 Calories.  Interestingly, I have ridden more hours in a month.  In June of 2016 I cycled 49 hours, 13 minutes, but only logged 469.3 miles.  Clearly my average speed has increased.  In addition I have noticed my average heart rate has dropped and, naturally, so has the number of calories burned.  I'm doing less work for more miles, thus becoming more efficient.  It is interesting to see how quickly my heart rate drops while waiting for a stop light or crossing a busy road.
   The only reason I was able to ride this many miles in October was the unseasonably warm weather we have been having.  The last week has been more typical and much harder to ride and it will only get cooler.  So I expect this new person record of mine will stand for at least a few months.  It will be hard to beat, but not impossible.  I'm not exhausted like I was last month.  Part of this maybe allowing myself to eat whatever I want, and trying to put a lot of protein into that diet.  My weight has been fairly stable and my fatigue has weened.  When the spring arrives I might just have to push myself again to see what I can get my body to do. 
Cat Burgler

Cat Burgler

   Biked downtown to have breakfast/lunch/dinner.  Although it was forecast to reach 50°F (10°C) this afternoon the temperature only hit 45°F (7°C).  Even with the wind however it was warm enough to just ride in a long sleeve shirt.  I plan to ride tomorrow in similar temperatures but the wind is forecast to be higher.  I keep a sweatshirt strapped to my bike but I shouldn't need to use it.