In this post, we are going to test different algorithms in C++ for fast median calculation on a microcontroller like Arduino. The interest we have in calculating the median is that, in certain applications, it is a better statistic compared to others such as the mean.
In the multiple sampling post, we saw that to improve the accuracy in data acquisition, the usual practice is to take multiple measurements and combine the different elements in some calculation.
The first candidate that surely comes to mind is to use the mean, as we will see when doing a moving average filter. However, the mean has the problem of being affected by noise and especially by spurious points.
Imagine that we take ten measurements from a sensor, and we get 9 values between 19.7 and 19.9. But one of the measurements gives an anomalous point with a value of 1000. A single value has caused the mean to go to almost 58, ruining the entire series.
If you want a more everyday and closer example, 10 people in a room earn an average of 1000 euros. But it may be that 9 are earning 100 euros and the last one 9100. That is the problem of the mean.
In these circumstances, the use of a statistic such as the median can be more representative of the trend. In fact, “robust statistics,” an alternative representation to classical statistics focused on the robustness of the data, uses the median instead of the mean as an indicator of trend.
One of the reasons why the median has traditionally been overshadowed by the mean as an indicator is not its greater statistical significance but, simply, the latter is simple and quick to calculate. Furthermore, it is a function that we can manipulate mathematically and derive to use in optimization solutions.
In contrast, calculating the median is a costly process. Additionally, it is not a function with which we can operate, and it is not derivable. However, the median has the advantage of being more robust to the appearance of noise, especially from spurious points.
This makes the median very common in acquisition filters, where it shows better behavior especially in salt and pepper noise filtering. We will see this when applying median filtering.
It is commonly considered that calculating the median requires sorting the elements of the sample and taking the middle term. However, there are more efficient algorithms for median calculation.
Here we are going to test three different algorithms with two samples of 100 and 500 elements on an Arduino Nano 16Mhz. As we already mentioned in the QuickSort post, due to its limitations, it is not the most suitable machine for calculations, and if we need these calculations, we should consider a more powerful processor. But doing the tests with 10 elements would be boring!
QuickSort Algorithm
The base case study is to sort the input vector and take the middle point. For this, we will use the QuickSort algorithm that we saw in the last post as one of the fastest algorithms for sorting a vector.
int values100[] = { 3, 53, 70, 56, 18, 85, 27, 14, 37, 94, 9, 55, 40, 60, 52, 61, 15, 65, 13, 8, 57, 97, 69, 4, 35, 82, 22, 73, 59, 68, 78, 24, 21, 36, 71, 80, 74, 39, 17, 12, 29, 76, 49, 51, 30, 90, 88, 2, 84, 50, 62, 28, 77, 43, 5, 16, 58, 26, 32, 34, 1, 75, 66, 95, 38, 89, 67, 87, 100, 54, 92, 81, 25, 83, 46, 33, 23, 45, 96, 99, 79, 48, 11, 31, 7, 6, 19, 91, 93, 44, 47, 98, 86, 41, 63, 20, 72, 10, 42, 64 };
int values100Length = sizeof(values100) / sizeof(int);
int values500[] = { 404, 267, 446, 214, 45, 149, 475, 496, 233, 78, 248, 307, 95, 431, 479, 445, 181, 370, 458, 476, 371, 122, 231, 74, 8, 392, 355, 397, 426, 125, 15, 159, 172, 369, 441, 318, 203, 399, 249, 225, 457, 351, 462, 184, 384, 100, 265, 244, 32, 499, 448, 29, 412, 447, 110, 473, 12, 414, 311, 301, 56, 84, 243, 378, 210, 217, 165, 10, 79, 374, 337, 52, 373, 395, 30, 126, 116, 280, 313, 474, 157, 6, 467, 459, 381, 129, 482, 13, 179, 167, 72, 68, 112, 194, 205, 97, 342, 142, 4, 418, 22, 440, 430, 364, 82, 483, 158, 198, 124, 259, 20, 312, 241, 254, 456, 361, 5, 245, 281, 376, 461, 274, 219, 348, 235, 23, 328, 2, 136, 291, 455, 302, 107, 415, 393, 43, 427, 211, 223, 168, 340, 87, 286, 133, 228, 354, 182, 204, 67, 419, 63, 270, 463, 60, 49, 358, 362, 102, 330, 242, 406, 108, 221, 83, 300, 363, 166, 290, 389, 436, 263, 34, 487, 377, 106, 491, 434, 257, 207, 417, 47, 379, 343, 500, 339, 403, 390, 61, 495, 262, 128, 132, 293, 94, 69, 143, 279, 375, 421, 109, 237, 310, 432, 218, 161, 150, 470, 200, 121, 464, 494, 443, 466, 252, 33, 105, 173, 344, 275, 388, 289, 333, 409, 452, 118, 315, 489, 283, 433, 442, 439, 114, 334, 229, 304, 175, 253, 216, 236, 256, 70, 169, 321, 365, 405, 366, 91, 380, 37, 212, 429, 336, 141, 308, 90, 492, 31, 460, 324, 387, 156, 120, 24, 183, 401, 81, 51, 288, 3, 367, 246, 498, 39, 386, 36, 192, 352, 292, 451, 294, 50, 326, 345, 76, 319, 360, 335, 306, 48, 239, 309, 468, 331, 226, 385, 347, 295, 44, 89, 497, 438, 332, 297, 346, 25, 199, 485, 469, 55, 402, 193, 284, 264, 135, 7, 14, 53, 197, 88, 201, 232, 258, 234, 481, 255, 9, 186, 59, 162, 18, 98, 93, 154, 73, 327, 278, 230, 131, 145, 26, 465, 103, 220, 19, 316, 177, 407, 146, 42, 153, 144, 99, 117, 28, 62, 148, 282, 453, 137, 208, 424, 450, 1, 477, 164, 368, 119, 21, 396, 127, 484, 277, 130, 437, 111, 206, 64, 391, 322, 151, 372, 188, 191, 57, 160, 383, 411, 180, 104, 16, 147, 170, 285, 350, 394, 71, 190, 54, 251, 261, 272, 320, 196, 338, 425, 227, 178, 420, 123, 174, 276, 480, 138, 80, 486, 454, 490, 435, 400, 77, 444, 323, 140, 171, 139, 195, 260, 314, 92, 17, 449, 222, 187, 296, 96, 40, 58, 423, 152, 11, 269, 359, 329, 422, 410, 155, 250, 101, 240, 213, 478, 273, 189, 27, 471, 356, 416, 238, 35, 41, 398, 113, 268, 46, 215, 85, 488, 325, 163, 202, 247, 341, 382, 299, 185, 176, 224, 472, 115, 349, 271, 303, 287, 408, 428, 65, 134, 75, 305, 66, 298, 357, 38, 266, 353, 209, 493, 317, 86, 413 };
int values500Length = sizeof(values500) / sizeof(int);
void setup()
{
Serial.begin(115200);
Serial.println("Median 100 values");
long timeCount = micros();
QuickSortAsc(values100, 0, values100Length - 1);
timeCount = micros() - timeCount;
Serial.println(values100[values100Length/2-1]);
Serial.println();
Serial.print(timeCount);
Serial.println("us");
Serial.println("");
Serial.println("Median 500 values");
timeCount = micros();
QuickSortAsc(values500, 0, values500Length - 1);
timeCount = micros() - timeCount;
Serial.println(values500[values500Length/2-1]);
Serial.println();
Serial.print(timeCount);
Serial.println("us");
}
void loop()
{
}
void QuickSortAsc(int* arr, const int left, const int right)
{
int i = left, j = right;
int tmp;
int pivot = arr[(left + right) / 2];
while (i <= j)
{
while (arr[i]<pivot) i++;
while (arr[j]>pivot) j--;
if (i <= j)
{
tmp = arr[i];
arr[i] = arr[j];
arr[j] = tmp;
i++;
j--;
}
};
if (left<j)
QuickSortAsc(arr, left, j);
if (i<right)
QuickSortAsc(arr, i, right);
}These are the results.
Logically, the times are identical to those obtained in the QuickSort post, 1540us for the sample of 100 elements, and 9500us for the sample of 500 elements.
QuickSelect Algorithm
The QuickSelect algorithm allows obtaining the k-th element of a list. It is based on QuickSort developed by the same author (Tony Hoare) and shares the partition and pivot system. However, by not having to sort the entire list, the order of QuickSelect goes from O(n log n) of QuickSort to O(n) in QuickSelect, maintaining in the worst-case O(n2).
#define ELEM_SWAP(a,b) { register float t=(a);(a)=(b);(b)=t; }
int values100[] = { 3, 53, 70, 56, 18, 85, 27, 14, 37, 94, 9, 55, 40, 60, 52, 61, 15, 65, 13, 8, 57, 97, 69, 4, 35, 82, 22, 73, 59, 68, 78, 24, 21, 36, 71, 80, 74, 39, 17, 12, 29, 76, 49, 51, 30, 90, 88, 2, 84, 50, 62, 28, 77, 43, 5, 16, 58, 26, 32, 34, 1, 75, 66, 95, 38, 89, 67, 87, 100, 54, 92, 81, 25, 83, 46, 33, 23, 45, 96, 99, 79, 48, 11, 31, 7, 6, 19, 91, 93, 44, 47, 98, 86, 41, 63, 20, 72, 10, 42, 64 };
int values100Length = sizeof(values100) / sizeof(int);
int values500[] = { 404, 267, 446, 214, 45, 149, 475, 496, 233, 78, 248, 307, 95, 431, 479, 445, 181, 370, 458, 476, 371, 122, 231, 74, 8, 392, 355, 397, 426, 125, 15, 159, 172, 369, 441, 318, 203, 399, 249, 225, 457, 351, 462, 184, 384, 100, 265, 244, 32, 499, 448, 29, 412, 447, 110, 473, 12, 414, 311, 301, 56, 84, 243, 378, 210, 217, 165, 10, 79, 374, 337, 52, 373, 395, 30, 126, 116, 280, 313, 474, 157, 6, 467, 459, 381, 129, 482, 13, 179, 167, 72, 68, 112, 194, 205, 97, 342, 142, 4, 418, 22, 440, 430, 364, 82, 483, 158, 198, 124, 259, 20, 312, 241, 254, 456, 361, 5, 245, 281, 376, 461, 274, 219, 348, 235, 23, 328, 2, 136, 291, 455, 302, 107, 415, 393, 43, 427, 211, 223, 168, 340, 87, 286, 133, 228, 354, 182, 204, 67, 419, 63, 270, 463, 60, 49, 358, 362, 102, 330, 242, 406, 108, 221, 83, 300, 363, 166, 290, 389, 436, 263, 34, 487, 377, 106, 491, 434, 257, 207, 417, 47, 379, 343, 500, 339, 403, 390, 61, 495, 262, 128, 132, 293, 94, 69, 143, 279, 375, 421, 109, 237, 310, 432, 218, 161, 150, 470, 200, 121, 464, 494, 443, 466, 252, 33, 105, 173, 344, 275, 388, 289, 333, 409, 452, 118, 315, 489, 283, 433, 442, 439, 114, 334, 229, 304, 175, 253, 216, 236, 256, 70, 169, 321, 365, 405, 366, 91, 380, 37, 212, 429, 336, 141, 308, 90, 492, 31, 460, 324, 387, 156, 120, 24, 183, 401, 81, 51, 288, 3, 367, 246, 498, 39, 386, 36, 192, 352, 292, 451, 294, 50, 326, 345, 76, 319, 360, 335, 306, 48, 239, 309, 468, 331, 226, 385, 347, 295, 44, 89, 497, 438, 332, 297, 346, 25, 199, 485, 469, 55, 402, 193, 284, 264, 135, 7, 14, 53, 197, 88, 201, 232, 258, 234, 481, 255, 9, 186, 59, 162, 18, 98, 93, 154, 73, 327, 278, 230, 131, 145, 26, 465, 103, 220, 19, 316, 177, 407, 146, 42, 153, 144, 99, 117, 28, 62, 148, 282, 453, 137, 208, 424, 450, 1, 477, 164, 368, 119, 21, 396, 127, 484, 277, 130, 437, 111, 206, 64, 391, 322, 151, 372, 188, 191, 57, 160, 383, 411, 180, 104, 16, 147, 170, 285, 350, 394, 71, 190, 54, 251, 261, 272, 320, 196, 338, 425, 227, 178, 420, 123, 174, 276, 480, 138, 80, 486, 454, 490, 435, 400, 77, 444, 323, 140, 171, 139, 195, 260, 314, 92, 17, 449, 222, 187, 296, 96, 40, 58, 423, 152, 11, 269, 359, 329, 422, 410, 155, 250, 101, 240, 213, 478, 273, 189, 27, 471, 356, 416, 238, 35, 41, 398, 113, 268, 46, 215, 85, 488, 325, 163, 202, 247, 341, 382, 299, 185, 176, 224, 472, 115, 349, 271, 303, 287, 408, 428, 65, 134, 75, 305, 66, 298, 357, 38, 266, 353, 209, 493, 317, 86, 413 };
int values500Length = sizeof(values500) / sizeof(int);
void setup()
{
Serial.begin(115200);
Serial.println("Mediana 100 valores");
long timeCount = micros();
QuickSelectMedian(values100, values100Length);
timeCount = micros() - timeCount;
Serial.println(values100[values100Length/2-1]);
Serial.println();
Serial.print(timeCount);
Serial.println("us");
Serial.println("");
Serial.println("Mediana 500 valores");
timeCount = micros();
QuickSelectMedian(values500, values500Length);
timeCount = micros() - timeCount;
Serial.println(values500[values500Length/2-1]);
Serial.println();
Serial.print(timeCount);
Serial.println("us");
}
void loop()}
int QuickSelectMedian(int arr[], uint16_t n)
{
uint16_t low, high;
uint16_t median;
uint16_t middle, ll, hh;
low = 0; high = n - 1; median = (low + high) / 2;
for (;;)
{
if (high <= low)
return arr[median];
if (high == low + 1)
{
if (arr[low] > arr[high])
ELEM_SWAP(arr[low], arr[high]);
return arr[median];
}
middle = (low + high) / 2;
if (arr[middle] > arr[high])
ELEM_SWAP(arr[middle], arr[high]);
if (arr[low] > arr[high])
ELEM_SWAP(arr[low], arr[high]);
if (arr[middle] > arr[low])
ELEM_SWAP(arr[middle], arr[low]);
ELEM_SWAP(arr[middle], arr[low + 1]);
ll = low + 1;
hh = high;
for (;;)
{
do ll++; while (arr[low] > arr[ll]);
do hh--; while (arr[hh] > arr[low]);
if (hh < ll)
break;
ELEM_SWAP(arr[ll], arr[hh]);
}
ELEM_SWAP(arr[low], arr[hh]);
if (hh <= median)
low = ll;
if (hh >= median)
high = hh - 1;
}
return arr[median];
}The results are as follows.
The time for a sample of 100 elements using QuickSelect is 172us, 8.95 times faster than sorting the list. In the case of a sample of 500 elements, the time is 2676us, 3.55 times faster than performing list sorting.
Wirth Algorithm
The algorithm developed by Niklaus Wirth is a variation of Hoare’s QuickSelect algorithm that allows for a simpler and more efficient implementation. It was introduced in 1976 in the book ‘Algorithms + data structures = programs’
int values100[] = { 3, 53, 70, 56, 18, 85, 27, 14, 37, 94, 9, 55, 40, 60, 52, 61, 15, 65, 13, 8, 57, 97, 69, 4, 35, 82, 22, 73, 59, 68, 78, 24, 21, 36, 71, 80, 74, 39, 17, 12, 29, 76, 49, 51, 30, 90, 88, 2, 84, 50, 62, 28, 77, 43, 5, 16, 58, 26, 32, 34, 1, 75, 66, 95, 38, 89, 67, 87, 100, 54, 92, 81, 25, 83, 46, 33, 23, 45, 96, 99, 79, 48, 11, 31, 7, 6, 19, 91, 93, 44, 47, 98, 86, 41, 63, 20, 72, 10, 42, 64 };
int values100Length = sizeof(values100) / sizeof(int);
int values500[] = { 404, 267, 446, 214, 45, 149, 475, 496, 233, 78, 248, 307, 95, 431, 479, 445, 181, 370, 458, 476, 371, 122, 231, 74, 8, 392, 355, 397, 426, 125, 15, 159, 172, 369, 441, 318, 203, 399, 249, 225, 457, 351, 462, 184, 384, 100, 265, 244, 32, 499, 448, 29, 412, 447, 110, 473, 12, 414, 311, 301, 56, 84, 243, 378, 210, 217, 165, 10, 79, 374, 337, 52, 373, 395, 30, 126, 116, 280, 313, 474, 157, 6, 467, 459, 381, 129, 482, 13, 179, 167, 72, 68, 112, 194, 205, 97, 342, 142, 4, 418, 22, 440, 430, 364, 82, 483, 158, 198, 124, 259, 20, 312, 241, 254, 456, 361, 5, 245, 281, 376, 461, 274, 219, 348, 235, 23, 328, 2, 136, 291, 455, 302, 107, 415, 393, 43, 427, 211, 223, 168, 340, 87, 286, 133, 228, 354, 182, 204, 67, 419, 63, 270, 463, 60, 49, 358, 362, 102, 330, 242, 406, 108, 221, 83, 300, 363, 166, 290, 389, 436, 263, 34, 487, 377, 106, 491, 434, 257, 207, 417, 47, 379, 343, 500, 339, 403, 390, 61, 495, 262, 128, 132, 293, 94, 69, 143, 279, 375, 421, 109, 237, 310, 432, 218, 161, 150, 470, 200, 121, 464, 494, 443, 466, 252, 33, 105, 173, 344, 275, 388, 289, 333, 409, 452, 118, 315, 489, 283, 433, 442, 439, 114, 334, 229, 304, 175, 253, 216, 236, 256, 70, 169, 321, 365, 405, 366, 91, 380, 37, 212, 429, 336, 141, 308, 90, 492, 31, 460, 324, 387, 156, 120, 24, 183, 401, 81, 51, 288, 3, 367, 246, 498, 39, 386, 36, 192, 352, 292, 451, 294, 50, 326, 345, 76, 319, 360, 335, 306, 48, 239, 309, 468, 331, 226, 385, 347, 295, 44, 89, 497, 438, 332, 297, 346, 25, 199, 485, 469, 55, 402, 193, 284, 264, 135, 7, 14, 53, 197, 88, 201, 232, 258, 234, 481, 255, 9, 186, 59, 162, 18, 98, 93, 154, 73, 327, 278, 230, 131, 145, 26, 465, 103, 220, 19, 316, 177, 407, 146, 42, 153, 144, 99, 117, 28, 62, 148, 282, 453, 137, 208, 424, 450, 1, 477, 164, 368, 119, 21, 396, 127, 484, 277, 130, 437, 111, 206, 64, 391, 322, 151, 372, 188, 191, 57, 160, 383, 411, 180, 104, 16, 147, 170, 285, 350, 394, 71, 190, 54, 251, 261, 272, 320, 196, 338, 425, 227, 178, 420, 123, 174, 276, 480, 138, 80, 486, 454, 490, 435, 400, 77, 444, 323, 140, 171, 139, 195, 260, 314, 92, 17, 449, 222, 187, 296, 96, 40, 58, 423, 152, 11, 269, 359, 329, 422, 410, 155, 250, 101, 240, 213, 478, 273, 189, 27, 471, 356, 416, 238, 35, 41, 398, 113, 268, 46, 215, 85, 488, 325, 163, 202, 247, 341, 382, 299, 185, 176, 224, 472, 115, 349, 271, 303, 287, 408, 428, 65, 134, 75, 305, 66, 298, 357, 38, 266, 353, 209, 493, 317, 86, 413 };
int values500Length = sizeof(values500) / sizeof(int);
#define ELEM_SWAP(a,b) { register float t=(a);(a)=(b);(b)=t; }
#define median(a,n) wirth_kth_smallest(a,n,(((n)&1)?((n)/2):(((n)/2)-1)))
void setup()
{
Serial.begin(115200);
int med;
Serial.println("Mediana 100 valores");
long timeCount = micros();
med = median(values100, values100Length);
timeCount = micros() - timeCount;
Serial.println(med);
Serial.println();
Serial.print(timeCount);
Serial.println("us");
Serial.println("");
Serial.println("Mediana 500 valores");
timeCount = micros();
med = median(values500, values500Length);
timeCount = micros() - timeCount;
Serial.println(med);
Serial.println();
Serial.print(timeCount);
Serial.println("us");
}
void loop()
{
}
int wirth_kth_smallest(int a[], int n, int k)
{
int i, j, l, m;
int x;
l = 0;
m = n - 1;
while (l<m)
{
x = a[k];
i = l;
j = m;
do {
while (a[i]<x) i++;
while (x<a[j]) j--;
if (i <= j)
{
ELEM_SWAP(a[i], a[j]);
i++;
j--;
}
} while (i <= j);
if (j<k) l = i;
if (k<i) m = j;
}
return a[k];
}The results are as follows.
The time for a sample of 100 elements is 172us, 8.95 times faster than QuickSort and exactly the same as the QuickSelect algorithm. In the case of a sample of 500 elements, it’s 1480us, 6.42 times faster than sorting the list, and 1.80 times faster than QuickSelect.
QuickMedian in a Library
What if we put it into a library to make it more convenient to use? Of course, here you have the algorithms in a QuickMedian library for Arduino. librería de QuickMedian para Arduino. Enjoy!
Descarga el código
Todo el código de esta entrada está disponible para su descarga en Github. 