mborgerson
Well-known member
During the thread on megasample collection with the T4X https://forum.pjrc.com/threads/60299-1-MSPS-on-a-T4-Is-this-possible we discussed various ways to evaluate noise and the best estimate of a stable signal on recorded ADC data. There were proponents of Pk-Pk noise, mean and standard deviation and median. I decided to evaluate the possibilities and results with both simulated normally-distributed noise and with real-world ADC data. I found that in most cases the mean and median were the same as a measure of central tendency for almost all the case with 12-bit data.
To accomplish this testing, I wrote the attached application that I have been running on a T4.1. The app has the following features:
* You can vary the ADC sampling and conversion speed.
* You can switch between ADC data and simulated data. With simulated data, what used to control ADC sampling sets the standard deviation of the simulated data.
* As part of the results there is a display of histogram data from the data set.
* There is a very fast median calculator---about 10X the speed of calculating the median by selecting the middle entry in quicksorted data.
* There is a very remedial plotting function that you can use with the Arduino Serial Plotter to plot the histogram data.
* You can add fake noise spikes to the data to see how they affect the statistics.
Here is the code:
To accomplish this testing, I wrote the attached application that I have been running on a T4.1. The app has the following features:
* You can vary the ADC sampling and conversion speed.
* You can switch between ADC data and simulated data. With simulated data, what used to control ADC sampling sets the standard deviation of the simulated data.
* As part of the results there is a display of histogram data from the data set.
* There is a very fast median calculator---about 10X the speed of calculating the median by selecting the middle entry in quicksorted data.
* There is a very remedial plotting function that you can use with the Arduino Serial Plotter to plot the histogram data.
* You can add fake noise spikes to the data to see how they affect the statistics.
Here is the code:
Code:
/******************************************************************
Program to demonstrate some sampling statistics to evaluate noise
in ADC collection.
Included is a fast algorithm for calculation of median and histogram
data for 12-bit integer ADC counts.
You can also add fake noise spikes to the data to see the effect
on the mean, median, std. deviation and histogram
M. Borgerson 8/18/2020
*********************************************************************/
#include <ADC.h>
const int ledpin = 13;
#define SIMMEAN 2200 // asumed mean for simulated data
#define SPIKEAMPLITUDE 20 // spike amplitude in counts--spikes always positive
#define SPIKEINTERVAL 100 // number of samples between spikes
const char compileTime [] = "T4.1 Noise Measurement Test Compiled on " __DATE__ " " __TIME__;
//NOTE: If you less than 65535 samples, you should change the bins array to uint16_t
// to save some space. Teensies have lots of memory, so 16KB for bins isn't an issue.
// However, if you collect 16-bit data, 65536 bins might be a problem.
// if you want a LOT of samples on a T4.1 with PSRAM increase this to 4 million and put
// idata in EXTMEM
#define MAXVALS 40000
#define MAXBINS 4096
#define LEDON digitalWriteFast(ledpin, HIGH);
#define LEDOFF digitalWriteFast(ledpin, LOW);
uint16_t idata[MAXVALS];
uint32_t bins[MAXBINS]; // sufficient for up to LOTS of 12-bit data
// instantiate a new ADC object
ADC *adc = new ADC(); // adc object;
bool simdata = true;
bool noisespikes = false;
void setup() {
pinMode(ledpin, OUTPUT);
pinMode(A0, INPUT_DISABLE); // disable digital keeper resistors
Serial.begin(9600);
delay(500); // wait for Serial to open
Serial.println();
Serial.println(compileTime);
Serial.printf("CPU Frequency: %lu\n", F_CPU_ACTUAL);
srand(micros()); // use micros() to set seed
}
void loop() {
char ch;
uint16_t cval;
if (Serial.available()) {
ch = Serial.read();
switch (ch) {
case 'b' : // for checking serial plotter input
LEDON
delay(1000);
LEDOFF
break;
case 's' :
Serial.println("Simulated data will be used.");
simdata = true;
break;
case 'a' :
Serial.println("Analog data will be used.");
simdata = false;
break;
case 'n' :
noisespikes = !noisespikes;
Serial.print("Noise spikes are " );
if (noisespikes) {
Serial.println(" ON");
} else {
Serial.println(" OFF");
}
break;
case '1' :
case '2' :
case '3' :
case '4' :
case '5' :
cval = ch & 0x07;
GenData(idata, MAXVALS, cval); // cval controls std for simulation, or ADC speed
HistoStats(idata, MAXVALS);
break;
case 'g' :
HistoGraph(bins);
LEDON
delay(100);
HistoGraph(bins);
LEDOFF
break;
} // end of switch(ch)
} // end of if(Serial.available....
} // end of loop()
// send bin data to Serial() so that Serial Minitor can plot it
// Doesn't plot any zero values, so low-count tails may be messed up
#define BINWIDTH 10
void HistoGraph(uint32_t gbins[]) {
uint16_t i, j;
// Serial.println("You have 10 seconds to switch to graph mode <Ctl-shift-L>)");
delay(10000);
for (i = 0; i < MAXBINS; i++) {
if (gbins[i] > 0) {
for (j = 0 ; j < BINWIDTH; j++) {
Serial.println(gbins[i]);
delay(10);
}
} // end of if(gbins[i] > 0)
} // end of for(i=0...
}
// Do some stats on values in idat
void HistoStats(uint16_t idat[], uint32_t numvals) {
uint32_t startmicro, dmicro;
uint16_t median, imin, imax;
float mean, std;
Serial.println("\nComputing statistics.");
delay(5);
MeanEtc(idat, numvals, &mean, &std, &imin, &imax);
Serial.printf("\nmean : % 8.4f ", mean);
Serial.printf("std. dev: % 8.4f ", std);
Serial.printf("min = % u max = % u\n \n", imin, imax);
Serial.println();
startmicro = micros();
median = binMedian(idat, numvals);
dmicro = micros() - startmicro;
Serial.printf("%lu sample bin Median = % u ", numvals, median);
Serial.printf(" took % 4lu microseconds\n", dmicro);
startmicro = micros();
qsort(idat, numvals, sizeof(uint16_t), cmpfunc);
median = idat[numvals / 2];
dmicro = micros() - startmicro;
Serial.printf("%lu sample QSort Median = % u ", numvals, median);
Serial.printf(" took % 4lu microseconds\n\n", dmicro);
// Now show histogram from bins data
ShowHisto(bins, imin, imax);
}
void GenData(uint16_t dat[], uint32_t nvals, uint16_t istd) {
if (simdata) {
SimGenData(dat, nvals, SIMMEAN, istd);
} else {
ADCGenData(dat, nvals, istd); // mean set by ADC input, istd controls speed setting
}
}
// do the easy stats
void MeanEtc(uint16_t dat[], uint32_t nvals, float *mnptr, float *stptr, uint16_t *minp, uint16_t *maxp) {
uint32_t i;
uint16_t imin, imax, ival;
float var, fv, s1, s2, mean;
double sum, sumsq; // with 4 million samples, even 32-bit floats have rounding problems!
// NOTE SUMSQ gets LARGE: suppose the values are ~4000 and you sum 100,000 of them then the sumsq
// is about (4x10^3) * (4x10^3) * (1x10^5) or 1.6x10^12 which is way too big for a uint32_t
sum = 0.0;
sumsq = 0.0;
imin = 65535;
imax = 0;
// first get sum and calculate mean
for (i = 0; i < nvals; i++) {
ival = dat[i];
if (ival < imin) imin = ival;
if (ival > imax) imax = ival;
sum += (float)ival;
sumsq += (float)ival * (float)ival;
}
mean = sum / nvals;
// now rcalculate variance in tw-pass algorithm to minimize roundoff error
s1 = 0; s2 = 0;
for (i = 0; i < nvals; i++) {
fv = (float)dat[i] - mean;
s1 += (fv * fv );
s2 += fv;
}
var = (s1 - s2 * s2 / nvals) / (nvals - 1);
*mnptr = mean;
//Serial.printf("Variance: %8.3ef s1: %8.3e s2^2/n: %8.3e\n", var, s1, s2*s2/nvals);
*stptr = sqrt(var);
*minp = imin;
*maxp = imax;
}
// Since this is aimed at evaluating noise of 12-bit ADC data,
// the generated numbers are constrained to values from 0 to 4095
// The numbers are also constrained to mean +/- 4 *std dev
void SimGenData(uint16_t dat[], uint32_t nvals, uint16_t imean, uint16_t istd) {
float x, fmean, fstd;
uint16_t ival;
uint32_t i;
bool badresult;
fmean = imean;
fstd = istd;
if (nvals > MAXVALS) nvals = MAXVALS; // stay within targt array
Serial.println("\n\nGenerating Random Normally distributed data.");
for (i = 0; i < nvals; i++) {
LEDON
do {
badresult = false;
x = fastrnorm(fmean, fstd);
ival = floor(x);
// if noisespikes are on add a spike every 100 samples
if (ival > MAXBINS - 1) badresult = true;
if (ival > (fmean + 6 * fstd)) badresult = true;
if (ival < (fmean - 6 * fstd)) badresult = true;
} while (badresult); // fall through when result in bounds
if (noisespikes) {
if ((i % SPIKEINTERVAL) == (SPIKEINTERVAL - 1)) {
ival += SPIKEAMPLITUDE;
}
}
if (ival > MAXBINS - 1) ival = MAXBINS - 1; // even spikes can get clipped!
LEDOFF
dat[i] = ival;
} // end of for(i=0...
}
// Collect ADC Data from A0. dmean is ignored. dstd
// controls ADC speed settings. No fancy timing control--
// just collect analog data as fast as analogRead will
// provide it.
void ADCGenData(uint16_t dat[], uint32_t nvals, uint16_t spd) {
uint32_t i, startmicro, dmicro;
uint16_t ival;
startmicro = micros();
if (nvals > MAXVALS) nvals = MAXVALS;
Serial.println("\n\nReading ADC data.");
SetADCSpeed(spd);
LEDON
for (i = 0; i < nvals; i++) {
// straightforward blocking analog read.
ival = adc->adc0->analogRead(A0);
// if noisespikes are on add a spike every 100 samples
if (noisespikes) {
if ((i % SPIKEINTERVAL) == (SPIKEINTERVAL - 1)) {
ival += SPIKEAMPLITUDE;
}
}
if (ival > 4095) ival = 4095;
dat[i] = ival;
} // end of for(i=0...
LEDOFF
dmicro = micros() - startmicro;
Serial.printf("Collecting % lu samples took % 4lu microseconds\n", nvals, dmicro);
}
// Set ADC Speed from low to very high
void SetADCSpeed(uint16_t spd) {
switch (spd) {
case 1:
Serial.println("ADC Speed set to 12 - bit LOW_SPEED.");
adc->adc0->setResolution(12); // set bits of resolution
adc->adc0->setConversionSpeed(ADC_CONVERSION_SPEED::LOW_SPEED);
adc->adc0->setSamplingSpeed(ADC_SAMPLING_SPEED::LOW_SPEED);
break;
case 2:
Serial.println("ADC Speed set to 12 - bit MED_SPEED.");
adc->adc0->setResolution(12); // set bits of resolution
adc->adc0->setConversionSpeed(ADC_CONVERSION_SPEED::MED_SPEED);
adc->adc0->setSamplingSpeed(ADC_SAMPLING_SPEED::MED_SPEED);
break;
case 3:
Serial.println("ADC Speed set to 12 - bit HIGH_SPEED.");
adc->adc0->setResolution(12); // set bits of resolution
adc->adc0->setConversionSpeed(ADC_CONVERSION_SPEED::HIGH_SPEED);
adc->adc0->setSamplingSpeed(ADC_SAMPLING_SPEED::HIGH_SPEED);
break;
case 4:
Serial.println("ADC Speed set to 12 - bit VERY_HIGH range.");
adc->adc0->setResolution(12); // set bits of resolution
adc->adc0->setConversionSpeed(ADC_CONVERSION_SPEED::VERY_HIGH_SPEED);
adc->adc0->setSamplingSpeed(ADC_SAMPLING_SPEED::VERY_HIGH_SPEED);
break;
case 5:
Serial.println("ADC Speed set to 8 - bit HIGH_SPEED.");
adc->adc0->setResolution(8); // set bits of resolution
adc->adc0->setConversionSpeed(ADC_CONVERSION_SPEED::HIGH_SPEED);
adc->adc0->setSamplingSpeed(ADC_SAMPLING_SPEED::HIGH_SPEED);
break;
}
}
// finds median without recursion, comparison or sorting
// takes about 140mSec to find median of 100,000 records
uint16_t binMedian(uint16_t dat[], uint32_t numsamp) {
uint32_t i, bsum;
uint16_t val;
// empty the bins--they're on the stack and have random values
memset(bins, 0, sizeof(bins));
// Serial.printf("In binMedian, address of dat is % p\n", &dat);
// delay(5);
for (i = 0; i < numsamp; i++) {
val = dat[i];
if (val < MAXBINS) { // prevent writing past end of bins array
bins[val]++;
}
} // bins are filled
i = 0; bsum = 0;
do {
bsum += bins[i];
i++;
} while (bsum < numsamp / 2);
return i - 1;
}
void ShowHisto(uint32_t bins[], uint16_t imin, uint16_t imax) {
uint16_t i, lnum;
lnum = 0;
Serial.print("Histogram Data");
for (i = imin - 1; i < imax + 2; i++) {
if ((lnum % 10) == 0) Serial.printf("\n % 4d: ", i);
lnum++;
Serial.printf(" % 5lu ", bins[i]);
}
Serial.println();
}
// These are the data collection and simulation routines
inline int cmpfunc (const void * a, const void * b) {
uint16_t ia = *(uint16_t *)a;
uint16_t ib = *(uint16_t *)b;
return ( ia - ib);
}
// returns a random floating point value between 0.0 and 1.0
// we use the built-in rand() function for speed as we are
// not particularly picky about the characteristics of the number
float randf(void) {
float val;
val = (float)rand() / ((float)(RAND_MAX) - 0.5);
return val;
}
// convert floating rnd to normal curve
float fastrnorm(float mean, float stdDev)
{
static float spare;
static float u1;
static float u2;
static float s;
static bool isSpareReady = false;
if (isSpareReady)
{
isSpareReady = false;
return ((spare * stdDev) + mean);
} else {
do {
u1 = (randf() * 2) - 1;
u2 = (randf() * 2) - 1;
s = (u1 * u1) + (u2 * u2);
} while (s >= 1.0);
s = sqrtf(-2.0 * logf(s) / s);
spare = u2 * s;
isSpareReady = true;
return (mean + (stdDev * u1 * s));
}
}