Go Down

Topic: exponentially distributed pseudorandom pulse generator (Read 453 times) previous topic - next topic

zike

This is my very first microcontroller project; constructive criticism appreciated! Usable as is, but from here I will take a different tack, so I'll post this as final before striking out in a new direction. Big picture, the object is to simulate output of a radiation detector, such as a Geiger-Muller tube or scintillator. This supports test and calibration on a longer-term counter/discriminator project.

I used a $4 Nano clone and junk box parts to build one, it works as intended. If I get a chance I'll scan the schematic and append it. The method and performance are seriously torqued by limitations of the AVR chip; this post is evidently over length (9000 characters???), so I'll maybe add a comment later to discuss methods and constraints (if others are interested).

After a failed dalliance with an Arduino 101, I ordered a 32-bit 180 MHz Teensy 3.6 (with floating point unit) to play with. Well over 100x faster, makes this whole task completely trivial. Not looking back...

Thanks for reading!

Zike

Code: [Select]

/*
   Random-interval pulse generator for testing radiation counter/scaler devices.
   MEZ  14 April 2017

   

*/

#include <avr/pgmspace.h>

// store lookup table in progam flash
const unsigned int samples[] PROGMEM = {
  827,   1037,    571,    165,    462,    716,    598,    737,
  213,    204,    321,    341,    349,    665,    266,    246,
  1395,    430,    800,    613,   1808,    114,    179,    333,
  676,     38,   1374,    198,    492,   1351,   2186,     28,
  971,    169,    641,    297,    274,   1372,      9,    174,
  174,   1563,    272,   2382,   2651,   1418,   1969,    185,

// ...CONTINUES MANY PAGES...MAKE YOUR OWN...

  26,    146,    432,    844,   3562,   2307,    348,    415,
  206,      7,    414,    564,    484,   1043,    116,    152,
  223,    243,    314,    682,     20,    512,    554,     69,
  857,    323,    241,    116,   1016,    125,    950,    144
}; // 8192 values, mean = 599.26; MATLAB function exprnd(600,8192);

unsigned long adjustf         = 16 ;       // cal offset for execution time (fixed interval)
unsigned long adjust          = 60 ;       // cal offset for execution time (random interval)
unsigned long rseed           = 3 ;        // random number seed
unsigned long pulseWidth      = 10 ;       // microseconds
unsigned long tauMean         = 600 ;      // microseconds
unsigned long powerOffTime    = 3600 * 1000L; // time to die (ms)

long          multiplier      = 1 ;        //
long          tau             = tauMean;   // negative values must be permitted
int           buffySize       = 1024;

byte          fixedRate       = 0;         // mode
byte          pulseOutPin     = 7;
byte          pulseBarOutPin  = 6;
byte          batteryCheckPin = 3;         // resistive divider connected to Vin
byte          controlRatePin  = 1;         // potentiometer wiper
byte          ledPin          = 11;
byte          interruptPin    = 2;
byte          killSwitchPin   = 9;

//------------------------------B76543210
byte          outToggle       = B00000000; // for setting bit mask


//unsigned long cum             = 0;        // diagnostics only
//unsigned long tic             = 0;        // diagnostics only
//unsigned long toc             = 0;        // diagnostics only
//unsigned long toe             = 0;        // diagnostics only
//unsigned long total           = 0;        // diagnostics only
//int           b               = 0;        // diagnostics only


void setup() {

  pinMode(killSwitchPin,      OUTPUT);
  digitalWrite(killSwitchPin, HIGH);      // maintain board power after boot
  pinMode(interruptPin,       INPUT);
  attachInterrupt(digitalPinToInterrupt(interruptPin), notAsimovThird, RISING); // next button press triggers kill

  pinMode(ledPin,             OUTPUT);
  pinMode(pulseOutPin,        OUTPUT);
  pinMode(pulseBarOutPin,     OUTPUT);

  digitalWrite(ledPin,        HIGH);
  delay(100);
//  digitalWrite(ledPin,        LOW);
  digitalWrite(pulseOutPin,   LOW);
  digitalWrite(pulseBarOutPin, HIGH);
  buffySize = sizeof(samples) / 2;      // 2 bytes per
  rseed = rseed + (analogRead(A4));     // use open ADC noise to scramble PRNG seed (unnecessary, but free!)


  //  Serial.begin(115200);
  //  while (!Serial && millis() < 200);
  //  Serial.println("Random-Interval Pulse Generator v. 1.0.5"); Serial.println(" MEZ, 14 April 2017 ");

  checkBattery();

  multiplier = readPot(controlRatePin);
  if (multiplier < 0) {
    fixedRate = 1;
    multiplier = -multiplier;
  }

  digitalWrite(ledPin,        HIGH);
  bitSet(outToggle, pulseOutPin);         // output  toggle mask set
  bitSet(outToggle, pulseBarOutPin);      // !output toggle mask set

  //  for (int k = 0 ; k < buffySize ; k++) { // check normalization of lookkup table
  //    total = total + pgm_read_word_near(samples + k);
  //  }
  //  Serial.println(); Serial.print("sample mean = "); Serial.println(total/buffySize);
  //  Serial.print("which should be close to "); Serial.println(tauMean);
  //
  //  tic = micros();

} // end of setup()


void loop() {
  if (millis() > powerOffTime) {       // shut down to save battery
    notAsimovThird();
  }

  if (fixedRate != 0) {
    tau = (multiplier * tauMean - adjustf - pulseWidth);
  }

  else {
    unsigned long j = (zrand() % buffySize);  // pick a random index in [0, buffySize)
    unsigned long value = (unsigned long) pgm_read_word_near(samples + j); // retrieve from prog flash
    tau = (multiplier * value - adjust - pulseWidth);
  }

  tau = max(tau, 0);            // avoid badness (floor)
  holdOff(tau);
  firePulse();

  //  diagnostics:
  //    toc = micros() - tic;
  //    if (b < 1000) {
  //      b++;
  //      cum = cum + toc ;
  //    }
  //    if (b == 1000) {
  //      Serial.println(cum);
  //      b = 0; cum = 0;
  //    }
  //    tic = micros();

} // end of loop()


/*
    functions
*/

void holdOff(unsigned long bau) {
  //  delayMicroseconds() fails beyond 16383 us; split long delays into ms + remainder us
  if (bau > 9999 ) {
    delay(bau / 1000) ;
    delayMicroseconds(bau % 1000 ) ;
  }
  else {
    delayMicroseconds(bau);
  }
} // end of holdOff

void firePulse() {
  // digitalWrite() takes 5 us per bit; direct register xor < 1 us
  // *** Atmega328 only *** change register name for other uP

  PORTD = PORTD ^ outToggle ;             // start pulse
  delayMicroseconds(pulseWidth);
  PORTD = PORTD ^ outToggle ;             // stop pulse
 
} // end of firePulse

unsigned long zrand() {
  // basic PRNG in style of C "rand()" function (Kernighan & Ritchie, 2nd ed., p.43)

  rseed = rseed * 1103515245 + 12345;
  return (unsigned long)(rseed / 65536) % 32768;
} // end of zrand()

void flashStatus(int nFlashes) {
  // show parameters by blinking LED code

  for (int k = 0 ; k < nFlashes ;  k++) {
    digitalWrite(ledPin, LOW);
    delay(150);
    digitalWrite(ledPin, HIGH);
    delay(150);
  }
  delay(500);
} // end of flashStatus()

long readPot(int pin) {
  //   read potentiometer voltage, set parameter values, signal selections

  int control = analogRead(pin) >> 7; // take only 3  MSB (0-7) (10-bit ADC boards only!)
  switch (control) {
    case 0:
      { flashStatus(1); flashStatus(2); return 1000;  break; // 100 cpm  (1.67 Hz)   random
      }
    case 1:
      { flashStatus(1); flashStatus(3); return 100;   break; // 1k cpm   (16.7 Hz)   random
      }
    case 2:
      { flashStatus(1); flashStatus(4); return 10;    break;  // 10k cpm  (167 Hz)   random
      }
    case 3:
      { flashStatus(1); flashStatus(5); return 1;     break;  // 100k cpm (1.67 kHz) random
      }
    case 4:
      { flashStatus(2); flashStatus(5); return -1;    break;  // 100k cpm (1.67 kHz) fixed
      }
    case 5:
      { flashStatus(2); flashStatus(4); return -10;   break;  // 10k cpm  (167 Hz)   fixed
      }
    case 6:
      { flashStatus(2); flashStatus(3); return -100;  break;  // 1k cpm   (16.7 Hz)  fixed
      }
    case 7:
      { flashStatus(2); flashStatus(2); return -1000; break;  // 100 cpm  (1.67 Hz)  fixed
      }
    default :
      { flashStatus(7); return 1;
      }
  }
} // end of readPot()

void checkBattery() {
  // ((7.5V-.25V)/5.0V ) * (100k/(348k+100k)) = 331/1023 ; .25V is schottky diode fwd drop
  if (analogRead(batteryCheckPin) < 331) {
    flashStatus(9);
  }
} // end of checkBattery()


void notAsimovThird() {
  digitalWrite(killSwitchPin, LOW);           // goodbye cruel world
} // end of everything



jremington

I can't imagine why a program to generate a pseudorandom exponential distribution needs to be more than a few lines long.

What is the lookup table for?

zike

I can't imagine why a program to generate a pseudorandom exponential distribution needs to be more than a few lines long.

What is the lookup table for?
Yep, that's what I thought! Long story...

Number of pulses in a given interval should be Poisson-distributed; equivalently, delays between adjacent pulses should be distributed exponentially. Simple enough: after each pulse, wait a random number of microseconds, selected from an exponential distribution, before issuing the next.

I set a design goal of simulating source activity of 100, 1k, 10k and 100k counts per minute (CPM) within ± 5%.  At 100k CPM the average inter-pulse interval should be 600 us. In practical terms, that means there shouldn't be any bias (e.g., missed or delayed pulses) for intervals as short as 30 us or so.

Did not take long to discover, a 16 MHz 8-bit Arduino ain't up to it. For example, digitalWrite() takes 5 us; the random() call takes 41 us; and floating point natural log(), required to remap uniform randoms into exponentially distributed, takes over 200 us.

Fallback plan: pre-load a sample of pseudorandom delays generated offline with "representative" statistics (the Table). However any finite set repeated sequentially fails miserably as a test signal (perfectly periodic). Furthermore, the Atmega 328P only has 2k of memory; at best you could cram in 512 samples before running out of room.

The latter issue is partly solved by using PROGMEM, so the array of constants occupies program memory rather than RAM. This allowed about 8192 delay samples, enough for a reasonably convincing histogram.

The former problem is addressed by selecting random array indices each time, scrambling the choices so no pattern repeats. Even this is too slow with library random(); so I included a 2-line "standard" PRNG from Kernighan and Ritchie (about 2 us per call on this hardware).

Finally, I toggled the digital output bit registers directly, both to reduce execution delay and to make sure complementary output pulses (+ and -) are aligned temporally.

In the code quoted I cut out most of the pseudorandom array declaration for readability; you can generate your own with a simple script (or call exprand() in Matlab  or random.expovariate() in Python). These samples were generated with a mean of 600 us, for the slower rates I just multiply by a constant.

Shorter story...

This is so damn complicated because a 16 MHz AVR is helpless. The task becomes trivial with a suitable processor.  In the end it works, I checked all the modes with an HP counter and they're within about 3% if you integrate long enough, and scatter plots and histograms look legit.

But with the Teensy I can calculate the delays on the fly and even get a decade higher rate (i.e., 1M CPS). With only a few lines of code.

BTW after finishing this I found a paper describing virtually the same thing implemented on an FPGA:
http://dx.doi.org/10.4236/wjnst.2013.34019

Attaching scan of my lab book with the schematic. Nothing special. Maybe one fun thing, power on/off/auto timeout circuit using 2 diodes and a P mosfet; inspired by Ken Moffatt via All About Circuits forum.

jremington

Yep, that application would not be a good fit for a 16-20 MHz, 8 bit processor!

zike

Some pics just for a laugh. My cnc panel engraver is down for maintenance...  :smiley-lol:

MarkT

There's probably some way to coax that performance from a ATmega 328, with a combination
of direct port access, LSFR RNG and interrupt driven sampling.  Set up timer2 to interrupt every
20us or similar, tune the threshold for each particular rate you want, and test every sample period,
no exponential needed:

Code: [Select]

ISR (TIMER2_OVF_vect)
{
    PORTD &= ~MASK ;
    unsigned long rand = lsfr_rng() ;  // use fast shift/xor implementation
    if (rand < threshold)
        PORTD |= MASK ;
}
[ I will NOT respond to personal messages, I WILL delete them, use the forum please ]

zike

Yes that's an elegant method!  Effectively your ensemble of rand() calls is a direct analogy for the ensemble of unstable nuclei, each with a fixed decay probability per unit time.  The disadvantage is that on average you need  to sample a substantial number of random numbers per unit (rate)^-1.  (other way to say it; to get within 5% or so you need to resolve times as short as 1/20 the mean interval, so you have to sample 20x per mean interval).

So your LSFR_RNG() must be really fast! Is there already an implementation for the AVR? 

Thanks,

Like

ard_newbie

#7
Apr 27, 2017, 07:06 pm Last Edit: Apr 27, 2017, 07:10 pm by ard_newbie
Just for my information, is an exponentially distributed pseudorandom something like:

X = - log(1 - Y)/ K   with Y a random value ( 0<Y<1) and K a constant ?

And BTW, Arduino DUE has a True Random Generator controller TRNG which produces a 4 bytes word every 1 us.

jremington

Quote
X = - log(1 - Y)/ K   with Y a random value ( 0<Y<1) and K a constant
Yes, if Y is uniformly distributed.

zike

That's interesting, can you point me to documentation? I considered testing a DUE but they don't seem to be available now.

The teensy also evidently has a hardware-based "true" RNG, but not very well documented (yet).

Of course for the subject application, pseudorandom is totally fine (so long as it's efficient).

And BTW, Arduino DUE has a True Random Generator controller TRNG which produces a 4 bytes word every 1 us.


ard_newbie

#11
Apr 30, 2017, 05:43 am Last Edit: Apr 30, 2017, 05:45 am by ard_newbie
You can find TRNG controller documentation chap. 42 page 1309 of sam3x datasheet.

See below a snippet to produce a True Random 32-bit number every 1 us:


Code: [Select]

/************************************************************************************/
/*     The True Random Number Generator (TRNG) on the DUE passes the American NIST  */
/*     Special Publication 800-22 and Diehard Random Tests Suites.                  */
/************************************************************************************/
#define TRNG_Key (0x524E47)

volatile uint32_t Random;

void setup() {

  PMC->PMC_PCER1 |= PMC_PCER1_PID41;                        // TRNG controller power ON
  TRNG->TRNG_CR = TRNG_CR_ENABLE | TRNG_CR_KEY(TRNG_Key);
  TRNG->TRNG_IER = TRNG_IER_DATRDY;
  NVIC_SetPriority(TRNG_IRQn, 0xFF);                        // Give TRNG interrupt the lowest urgency (if necessary )
  NVIC_EnableIRQ(TRNG_IRQn);

}


void loop() {

// Do whatever you want here with Random

}

/*     Fires every 1 us   ****/
void TRNG_Handler() {

  TRNG->TRNG_ISR;             // Clear status register
  Random = TRNG->TRNG_ODATA;  // Get random value

}





And it gives you a 1 us timer !!

MarkT

Yes that's an elegant method!  Effectively your ensemble of rand() calls is a direct analogy for the ensemble of unstable nuclei, each with a fixed decay probability per unit time.  The disadvantage is that on average you need  to sample a substantial number of random numbers per unit (rate)^-1.  (other way to say it; to get within 5% or so you need to resolve times as short as 1/20 the mean interval, so you have to sample 20x per mean interval).

So your LSFR_RNG() must be really fast! Is there already an implementation for the AVR? 

Thanks,

Like
<<, >>, ^
[ I will NOT respond to personal messages, I WILL delete them, use the forum please ]

Go Up