From f8af760ba077a6e7b0be00aac17e2490e14d6d31 Mon Sep 17 00:00:00 2001 From: Jeff Epler Date: Fri, 6 Mar 2020 11:50:59 -0600 Subject: [PATCH] Add examples for upcoming ulab guide --- ulab_Crunch_Numbers_Fast/benchmark.py | 40 +++++ ulab_Crunch_Numbers_Fast/cluebarometer.py | 175 ++++++++++++++++++++++ ulab_Crunch_Numbers_Fast/cluepulse.py | 148 ++++++++++++++++++ ulab_Crunch_Numbers_Fast/ledwave.py | 76 ++++++++++ ulab_Crunch_Numbers_Fast/waterfall.py | 96 ++++++++++++ 5 files changed, 535 insertions(+) create mode 100644 ulab_Crunch_Numbers_Fast/benchmark.py create mode 100644 ulab_Crunch_Numbers_Fast/cluebarometer.py create mode 100644 ulab_Crunch_Numbers_Fast/cluepulse.py create mode 100644 ulab_Crunch_Numbers_Fast/ledwave.py create mode 100644 ulab_Crunch_Numbers_Fast/waterfall.py diff --git a/ulab_Crunch_Numbers_Fast/benchmark.py b/ulab_Crunch_Numbers_Fast/benchmark.py new file mode 100644 index 000000000..25afb1c59 --- /dev/null +++ b/ulab_Crunch_Numbers_Fast/benchmark.py @@ -0,0 +1,40 @@ +import time +import math +import ulab +import ulab.numerical + +def mean(values): + return sum(values) / len(values) + +def normalized_rms(values): + minbuf = int(mean(values)) + samples_sum = sum( + float(sample - minbuf) * (sample - minbuf) + for sample in values + ) + + return math.sqrt(samples_sum / len(values)) + +def normalized_rms_ulab(values): + minbuf = ulab.numerical.mean(values) + values = values - minbuf + samples_sum = ulab.numerical.sum(values * values) + return math.sqrt(samples_sum / len(values)) + + +# Instead of using sensor data, we generate some data +# The amplitude is 5000 so the rms should be around 5000/1.414 = 3536 +nums_list = [int(8000 + math.sin(i) * 5000) for i in range(100)] +nums_array = ulab.array(nums_list) + +def timeit(s, f, n=100): + t0 = time.monotonic_ns() + for _ in range(n): + x = f() + t1 = time.monotonic_ns() + r = (t1 - t0) * 1e-6 / n + print("%-20s : %8.3fms [result=%f]" % (s, r, x)) + +print("Computing the RMS value of 1000 numbers") +timeit("traditional", lambda: normalized_rms(nums_list)) +timeit("ulab", lambda: normalized_rms_ulab(nums_array)) diff --git a/ulab_Crunch_Numbers_Fast/cluebarometer.py b/ulab_Crunch_Numbers_Fast/cluebarometer.py new file mode 100644 index 000000000..3096b9a6b --- /dev/null +++ b/ulab_Crunch_Numbers_Fast/cluebarometer.py @@ -0,0 +1,175 @@ +import time + +import adafruit_bmp280 +import board +import displayio +import ulab +import ulab.filter + +# Blank the screen. Scrolling text causes unwanted delays. +d = displayio.Group() +board.DISPLAY.show(d) + +# Sampling rate: 16Hz +# Cutoff frequency: 0.16Hz +# Transition bandwidth: 0.16Hz +# Window type: Hamming +# Filter has 311 coefficients +taps = ulab.array([ + -0.000050679794726066, -0.000041099278318167, -0.000031279920668665, + -0.000021183486597150, -0.000010770285292045, +0.000000000000000000, + +0.000011167446754809, +0.000022770999889941, +0.000034847558259864, + +0.000047431049079466, +0.000060551498686721, +0.000074234108254511, + +0.000088498343199344, +0.000103357045109305, +0.000118815575023601, + +0.000134870996840645, +0.000151511309510219, +0.000168714736477861, + +0.000186449080596567, +0.000204671152403140, +0.000223326279275218, + +0.000242347902542027, +0.000261657269119383, +0.000281163223679941, + +0.000300762106756334, +0.000320337763510928, +0.000339761667195315, + +0.000358893160569588, +0.000377579817760222, +0.000395657928211038, + +0.000412953103529159, +0.000429281007152519, +0.000444448205872873, + +0.000458253141344113, +0.000470487218795955, +0.000480936009263626, + +0.000489380560741255, +0.000495598812776238, +0.000499367108150093, + +0.000500461794444300, +0.000498660907473236, +0.000493745927786584, + +0.000485503600706003, +0.000473727809671115, +0.000458221492033063, + +0.000438798585855176, +0.000415285995764155, +0.000387525565446236, + +0.000355376044004699, +0.000318715033091691, +0.000277440901501588, + +0.000231474653767861, +0.000180761739242710, +0.000125273788160487, + +0.000065010261293197, +0.000000000000000000, -0.000069697336247377, + -0.000143989957415198, -0.000222752767634882, -0.000305826338672358, + -0.000393016043088374, -0.000484091357342654, -0.000578785344322494, + -0.000676794323931742, -0.000777777739462615, -0.000881358226495441, + -0.000987121890034750, -0.001094618794499868, -0.001203363670049808, + -0.001312836837542114, -0.001422485353209744, -0.001531724372895900, + -0.001639938734420840, -0.001746484755374530, -0.001850692242341569, + -0.001951866706278179, -0.002049291777482158, -0.002142231812333790, + -0.002229934682745978, -0.002311634738053158, -0.002386555927898205, + -0.002453915073551964, -0.002512925274028313, -0.002562799432345805, + -0.002602753886341418, -0.002632012127569287, -0.002649808591023194, + -0.002655392497711921, -0.002648031731496151, -0.002627016731069257, + -0.002591664377536210, -0.002541321857718479, -0.002475370483091317, + -0.002393229444145817, -0.002294359479963247, -0.002178266442894981, + -0.002044504738458277, -0.001892680620886388, -0.001722455325210333, + -0.001533548017297868, -0.001325738543930948, -0.001098869965763655, + -0.000852850856865069, -0.000587657355512251, -0.000303334951952833, + +0.000000000000000001, +0.000322159059450752, +0.000662880773589522, + +0.001021830060775982, +0.001398597909331569, +0.001792701398335994, + +0.002203584045179127, +0.002630616483032971, +0.003073097469789485, + +0.003530255228366684, +0.004001249116626688, +0.004485171623483914, + +0.004981050686118591, +0.005487852321559077, +0.006004483564265146, + +0.006529795699742466, +0.007062587782654920, +0.007601610426384373, + +0.008145569849526276, +0.008693132163411565, +0.009242927883419039, + +0.009793556645595150, +0.010343592108937170, +0.010891587022627668, + +0.011436078436539264, +0.011975593032464911, +0.012508652552774892, + +0.013033779302562583, +0.013549501700820601, +0.014054359855790191, + +0.014546911139352909, +0.015025735735186426, +0.015489442135386880, + +0.015936672560369614, +0.016366108277098043, +0.016776474791055797, + +0.017166546887869318, +0.017535153501103896, +0.017881182383493146, + +0.018203584559716979, +0.018501378539810983, +0.018773654273367416, + +0.019019576825867947, +0.019238389759765797, +0.019429418204303113, + +0.019592071599501125, +0.019725846101288819, +0.019830326636332028, + +0.019905188596781104, +0.019950199166862841, +0.019965218274992248, + +0.019950199166862841, +0.019905188596781104, +0.019830326636332028, + +0.019725846101288819, +0.019592071599501125, +0.019429418204303113, + +0.019238389759765800, +0.019019576825867947, +0.018773654273367420, + +0.018501378539810983, +0.018203584559716979, +0.017881182383493149, + +0.017535153501103892, +0.017166546887869318, +0.016776474791055797, + +0.016366108277098043, +0.015936672560369614, +0.015489442135386881, + +0.015025735735186426, +0.014546911139352912, +0.014054359855790193, + +0.013549501700820601, +0.013033779302562583, +0.012508652552774890, + +0.011975593032464912, +0.011436078436539264, +0.010891587022627668, + +0.010343592108937174, +0.009793556645595150, +0.009242927883419041, + +0.008693132163411567, +0.008145569849526276, +0.007601610426384373, + +0.007062587782654920, +0.006529795699742466, +0.006004483564265146, + +0.005487852321559078, +0.004981050686118592, +0.004485171623483914, + +0.004001249116626688, +0.003530255228366684, +0.003073097469789486, + +0.002630616483032971, +0.002203584045179127, +0.001792701398335993, + +0.001398597909331569, +0.001021830060775982, +0.000662880773589522, + +0.000322159059450752, +0.000000000000000001, -0.000303334951952833, + -0.000587657355512251, -0.000852850856865070, -0.001098869965763655, + -0.001325738543930948, -0.001533548017297868, -0.001722455325210333, + -0.001892680620886389, -0.002044504738458278, -0.002178266442894981, + -0.002294359479963247, -0.002393229444145818, -0.002475370483091317, + -0.002541321857718479, -0.002591664377536210, -0.002627016731069256, + -0.002648031731496151, -0.002655392497711923, -0.002649808591023195, + -0.002632012127569288, -0.002602753886341418, -0.002562799432345805, + -0.002512925274028314, -0.002453915073551965, -0.002386555927898205, + -0.002311634738053157, -0.002229934682745978, -0.002142231812333790, + -0.002049291777482159, -0.001951866706278179, -0.001850692242341569, + -0.001746484755374531, -0.001639938734420841, -0.001531724372895900, + -0.001422485353209745, -0.001312836837542114, -0.001203363670049807, + -0.001094618794499867, -0.000987121890034750, -0.000881358226495441, + -0.000777777739462615, -0.000676794323931742, -0.000578785344322494, + -0.000484091357342654, -0.000393016043088375, -0.000305826338672358, + -0.000222752767634882, -0.000143989957415198, -0.000069697336247377, + +0.000000000000000000, +0.000065010261293198, +0.000125273788160487, + +0.000180761739242710, +0.000231474653767861, +0.000277440901501588, + +0.000318715033091691, +0.000355376044004699, +0.000387525565446236, + +0.000415285995764155, +0.000438798585855176, +0.000458221492033064, + +0.000473727809671115, +0.000485503600706004, +0.000493745927786584, + +0.000498660907473236, +0.000500461794444300, +0.000499367108150093, + +0.000495598812776237, +0.000489380560741255, +0.000480936009263626, + +0.000470487218795956, +0.000458253141344114, +0.000444448205872873, + +0.000429281007152519, +0.000412953103529159, +0.000395657928211038, + +0.000377579817760223, +0.000358893160569588, +0.000339761667195315, + +0.000320337763510927, +0.000300762106756335, +0.000281163223679941, + +0.000261657269119383, +0.000242347902542027, +0.000223326279275218, + +0.000204671152403140, +0.000186449080596567, +0.000168714736477861, + +0.000151511309510219, +0.000134870996840645, +0.000118815575023601, + +0.000103357045109305, +0.000088498343199344, +0.000074234108254511, + +0.000060551498686721, +0.000047431049079466, +0.000034847558259864, + +0.000022770999889941, +0.000011167446754809, +0.000000000000000000, + -0.000010770285292045, -0.000021183486597150, -0.000031279920668665, + -0.000041099278318167, -0.000050679794726066, +]) + +# How often we are going to poll the sensor (If you change this, you need +# to change the filter above and the integration time below) +dt = 62500000 # 16Hz, 62.5ms + +# Wait until after deadline_ns has passed +def sleep_deadline(deadline_ns): + while time.monotonic_ns() < deadline_ns: + pass + + +# Initialize our sensor +i2c = board.I2C() +sensor = adafruit_bmp280.Adafruit_BMP280_I2C(i2c) +sensor.standby_period = adafruit_bmp280.STANDBY_TC_1000 +# Disable in-sensor filtering, because we want to show how it's done in +# CircuitPython +sensor.iir_filter = adafruit_bmp280.IIR_FILTER_DISABLE +sensor.overscan_pressure = adafruit_bmp280.OVERSCAN_X1 + +# And our data structures +# The most recent data samples, equal in number to the filter taps +data = ulab.zeros(len(taps)) +t0 = deadline = time.monotonic_ns() +n = 0 +# Take an initial reading to subtract off later, so that the graph in mu +# accentuates the small short term changes in pressure rather than the large +# DC offset of around 980 +offset = sensor.pressure + +while True: + deadline += dt + sleep_deadline(deadline) + # Move the trace near the origin so small differences can be seen in the mu + # plot window .. you wouldn't do this subtraction step if you are really + # interested in absolute barometric pressure. + value = sensor.pressure - offset + if n == 0: + # The first time, fill the filter with the initial value + data = data + value + else: + # Otherwise, add it as the next sample + ulab.numerical.roll(data, 1) + data[-1] = value + filtered = ulab.numerical.sum(data * taps) + # Actually print every 10th value. This prints about 1.6 values per + # second. You can print values more quickly by removing the 'if' and + # making the print unconditional, or change the frequency of prints up + # or down by changing the number '10'. + if n % 10 == 0: + print((filtered, value)) + n += 1 diff --git a/ulab_Crunch_Numbers_Fast/cluepulse.py b/ulab_Crunch_Numbers_Fast/cluepulse.py new file mode 100644 index 000000000..a1963d6f5 --- /dev/null +++ b/ulab_Crunch_Numbers_Fast/cluepulse.py @@ -0,0 +1,148 @@ +import time + +import adafruit_apds9960.apds9960 +import board +import digitalio +import ulab +import ulab.filter + +# Blank the screen. Scrolling text causes unwanted delays. +import displayio +d = displayio.Group() +board.DISPLAY.show(d) + +# Filter computed at https://fiiir.com/ +# Sampling rate: 8Hz +# Cutoff freqency: 0.5Hz +# Transition bandwidth 0.25Hz +# Window type: Regular +# Number of coefficients: 31 +# Manually trimmed to 16 coefficients +taps = ulab.array([ + +0.861745279666917052/2, + -0.134728583242092248, + -0.124472980501612152, + -0.108421190967457198, + -0.088015688587190874, + -0.065052714580474319, + -0.041490993500537393, + -0.019246940463156042, + -0.000000000000000005, + +0.014969842582454691, + +0.024894596100322432, + +0.029569415718397409, + +0.029338562862396955, + +0.025020274838643962, + +0.017781854357373172, + +0.008981905549472832, +]) + +# How much reflected light is required before pulse sensor activates +# These values are triggered when I bring my finger within a half inch. +# The sensor works when the finger is pressed lightly against the sensor. +PROXIMITY_THRESHOLD_HI = 225 +PROXIMITY_THRESHOLD_LO = 215 + +# These constants control how much the sensor amplifies received light +APDS9660_AGAIN_1X = 0 +APDS9660_AGAIN_4X = 1 +APDS9660_AGAIN_16X = 2 +APDS9660_AGAIN_64X = 3 + +# How often we are going to poll the sensor (If you change this, you need +# to change the filter above and the integration time below) +dt = 125000000 # 8Hz, 125ms + +# Wait until after deadline_ns has passed +def sleep_deadline(deadline_ns): + while time.monotonic_ns() < deadline_ns: + pass + +# Compute a high resolution crossing-time estimate for the sample, using a +# linear model +def estimated_cross_time(y0, y1, t0): + m = (y1 - y0) / dt + return t0 + round(-y1 / m) + +i2c = board.I2C() +sensor = adafruit_apds9960.apds9960.APDS9960(i2c) +white_leds = digitalio.DigitalInOut(board.WHITE_LEDS) +white_leds.switch_to_output(False) + +def main(): + sensor.enable_proximity = True + while True: + # Wait for user to put finger over sensor + while sensor.proximity() < PROXIMITY_THRESHOLD_HI: + time.sleep(.01) + + # After the finger is sensed, set up the color sensor + sensor.enable_color = True + # This sensor integration time is just a little bit shorter than 125ms, + # so we should always have a fresh value when we ask for it, without + # checking if a value is available. + sensor.integration_time = 220 + # In my testing, 64X gain saturated the sensor, so this is the biggest + # gain value that works properly. + sensor.color_gain = APDS9660_AGAIN_4X + white_leds.value = True + + # And our data structures + # The most recent data samples, equal in number to the filter taps + data = ulab.zeros(len(taps)) + # The filtered value on the previous iteration + old_value = 1 + # The times of the most recent pulses registered. Increasing this number + # makes the estimation more accurate, but at the expense of taking longer + # before a pulse number can be computed + pulse_times = [] + # The estimated heart rate based on the recent pulse times + rate = None + # the number of samples taken + n = 0 + + # Rather than sleeping for a fixed duration, we compute a deadline + # in nanoseconds and wait for the new deadline time to arrive. This + # helps the long term frequency of measurements better match the desired + # frequency. + t0 = deadline = time.monotonic_ns() + # As long as their finger is over the sensor, capture data + while sensor.proximity() >= PROXIMITY_THRESHOLD_LO: + deadline += dt + sleep_deadline(deadline) + value = sum(sensor.color_data) # Combination of all channels + ulab.numerical.roll(data, 1) + data[-1] = value + # Compute the new filtered variable by applying the filter to the + # recent data samples + filtered = ulab.numerical.sum(data * taps) + + # We gathered enough data to fill the filters, and + # the light value crossed the zero line in the positive direction + # Therefore we need to record a pulse + if n > len(taps) and old_value < 0 and filtered >= 0: + # This crossing time is estimated, but it increases the pulse + # estimate resolution quite a bit. If only the nearest 1/8s + # was used for pulse estimation, the smallest pulse increment + # that can be measured is 7.5bpm. + cross = estimated_cross_time(old_value, filtered, deadline) + # store this pulse time (in seconds since sensor-touch) + pulse_times.append((cross - t0) * 1e-9) + # and maybe delete an old pulse time + del pulse_times[:-10] + # And compute a rate based on the last recorded pulse times + if len(pulse_times) > 1: + rate = 60/(pulse_times[-1]-pulse_times[0])*(len(pulse_times)-1) + old_value = filtered + + # We gathered enough data to fill the filters, so report the light + # value and possibly the estimated pulse rate + if n > len(taps): + print((filtered, rate)) + n += 1 + + # Turn off the sensor and the LED and go back to the top for another run + sensor.enable_color = False + white_leds.value = False + print() +main() diff --git a/ulab_Crunch_Numbers_Fast/ledwave.py b/ulab_Crunch_Numbers_Fast/ledwave.py new file mode 100644 index 000000000..54e838656 --- /dev/null +++ b/ulab_Crunch_Numbers_Fast/ledwave.py @@ -0,0 +1,76 @@ +import random + +import board +import neopixel +from _pixelbuf import wheel +import ulab +import ulab.filter + +# Customize your neopixel configuration here... +pixel_pin = board.D5 +num_pixels = 96 +pixels = neopixel.NeoPixel(pixel_pin, num_pixels, brightness=0.1, + auto_write=False, pixel_order=neopixel.RGB) + +ddt = ulab.array([1.,-2.,1.]) +def step(u, um, f, n, dx, dt, c): + dt2 = dt*dt + C2 = (c*dt/dx)**2 + deriv = ulab.filter.convolve(u, ddt)[1:-1] * C2 + up = -um + u * 2 + deriv + f * dt2 + up[0] = 0 + up[n-1] = 0 + + return up + +def main(): + # This precomputes the color palette for maximum speed + # You could change it to compute the color palette of your choice + w = [wheel(i) for i in range(256)] + + # This sets up the initial wave as a smooth gradient + u = ulab.zeros(num_pixels) + um = ulab.zeros(num_pixels) + f = ulab.zeros(num_pixels) + + slope = ulab.linspace(0, 256, num=num_pixels) + th = 1 + + # the first time is always random (is that a contradiction?) + r = 0 + + while True: + + # Some of the time, add a random new wave to the mix + # increase .15 to add waves more often + # decrease it to add waves less often + if r < .01: + ii = random.randrange(1, num_pixels-1) + # increase 2 to make bigger waves + f[ii] = (random.random() - .5) * 2 + + # Here's where to change dx, dt, and c + # try .2, .02, 2 for relaxed + # try 1., .7, .2 for very busy / almost random + u, um = step(u, um, f, num_pixels, .1, .02, 1), u + + v = u * 200000 + slope + th + for i, vi in enumerate(v): + # Scale up by an empirical value, rotate by th, and look up the color + pixels[i] = w[round(vi) % 256] + + # Take away a portion of the energy of the waves so they don't get out + # of control + u = u * .99 + + # incrementing th causes the wheel to slowly cycle even if nothing else is happening + th = (th + .25) % 256 + pixels.show() + + # Clear out the old random value, if any + f[ii] = 0 + + # and get a new random value + r = random.random() + +main() diff --git a/ulab_Crunch_Numbers_Fast/waterfall.py b/ulab_Crunch_Numbers_Fast/waterfall.py new file mode 100644 index 000000000..b9641727d --- /dev/null +++ b/ulab_Crunch_Numbers_Fast/waterfall.py @@ -0,0 +1,96 @@ +"""Waterfall FFT demo adapted from +https://teaandtechtime.com/fft-circuitpython-library/ +to work with ulab on Adafruit CLUE""" + +import array + +import board +import audiobusio +import displayio +import ulab +import ulab.fft +import ulab.vector + +display = board.DISPLAY + +# Create a heatmap color palette +palette = displayio.Palette(52) +for i, pi in enumerate((0xff0000, 0xff0a00, 0xff1400, 0xff1e00, + 0xff2800, 0xff3200, 0xff3c00, 0xff4600, + 0xff5000, 0xff5a00, 0xff6400, 0xff6e00, + 0xff7800, 0xff8200, 0xff8c00, 0xff9600, + 0xffa000, 0xffaa00, 0xffb400, 0xffbe00, + 0xffc800, 0xffd200, 0xffdc00, 0xffe600, + 0xfff000, 0xfffa00, 0xfdff00, 0xd7ff00, + 0xb0ff00, 0x8aff00, 0x65ff00, 0x3eff00, + 0x17ff00, 0x00ff10, 0x00ff36, 0x00ff5c, + 0x00ff83, 0x00ffa8, 0x00ffd0, 0x00fff4, + 0x00a4ff, 0x0094ff, 0x0084ff, 0x0074ff, + 0x0064ff, 0x0054ff, 0x0044ff, 0x0032ff, + 0x0022ff, 0x0012ff, 0x0002ff, 0x0000ff)): + palette[51-i] = pi + +class RollingGraph(displayio.TileGrid): + def __init__(self, scale=2): + # Create a bitmap with heatmap colors + self.bitmap = displayio.Bitmap(display.width//scale, + display.height//scale, len(palette)) + super().__init__(self.bitmap, pixel_shader=palette) + + self.scroll_offset = 0 + + def show(self, data): + y = self.scroll_offset + bitmap = self.bitmap + + board.DISPLAY.auto_refresh = False + offset = max(0, (bitmap.width-len(data))//2) + for x in range(min(bitmap.width, len(data))): + bitmap[x+offset, y] = int(data[x]) + + board.DISPLAY.auto_refresh = True + + self.scroll_offset = (y + 1) % self.bitmap.height + +group = displayio.Group(scale=3) +graph = RollingGraph(3) +fft_size = 128 + +# Add the TileGrid to the Group +group.append(graph) + +# Add the Group to the Display +display.show(group) + +# instantiate board mic +mic = audiobusio.PDMIn(board.MICROPHONE_CLOCK, board.MICROPHONE_DATA, + sample_rate=16000, bit_depth=16) + +#use some extra sample to account for the mic startup +samples_bit = array.array('H', [0] * (fft_size+3)) + +# Main Loop +def main(): + max_all = 1 + + while True: + mic.record(samples_bit, len(samples_bit)) + samples = ulab.array(samples_bit[3:]) + spectrogram1 = ulab.fft.spectrum(samples) + # spectrum() is always nonnegative, but add a tiny value + # to change any zeros to nonzero numbers + spectrogram1 = ulab.vector.log(spectrogram1 + 1e-7) + spectrogram1 = spectrogram1[1:(fft_size//2)-1] + min_curr = min(spectrogram1) + max_curr = max(spectrogram1) + + if max_curr > max_all: + max_all = max_curr + else: + max_curr = max_curr-1 + + # Plot FFT + data = (spectrogram1 - min_curr) * (51. / (max_all - min_curr)) + graph.show(data) + +main()