diff --git a/ulab_Crunch_Numbers_Fast/benchmark.py b/ulab_Crunch_Numbers_Fast/benchmark.py
index fc1c421e7..384f23a5d 100644
--- a/ulab_Crunch_Numbers_Fast/benchmark.py
+++ b/ulab_Crunch_Numbers_Fast/benchmark.py
@@ -1,7 +1,6 @@
 import time
 import math
-import ulab
-import ulab.numerical
+import ulab.numpy
 
 def mean(values):
     return sum(values) / len(values)
@@ -17,15 +16,15 @@ def normalized_rms(values):
 
 def normalized_rms_ulab(values):
     # this function works with ndarrays only
-    minbuf = ulab.numerical.mean(values)
+    minbuf = ulab.numpy.mean(values)
     values = values - minbuf
-    samples_sum = ulab.numerical.sum(values * values)
+    samples_sum = ulab.numpy.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)
+nums_array = ulab.numpy.array(nums_list)
 
 def timeit(s, f, n=100):
     t0 = time.monotonic_ns()
@@ -38,5 +37,5 @@ def timeit(s, f, n=100):
 print("Computing the RMS value of 100 numbers")
 timeit("traditional", lambda: normalized_rms(nums_list))
 timeit("ulab, with ndarray, some implementation in python", lambda: normalized_rms_ulab(nums_array))
-timeit("ulab only, with list", lambda: ulab.numerical.std(nums_list))
-timeit("ulab only, with ndarray", lambda: ulab.numerical.std(nums_array))
+timeit("ulab only, with list", lambda: ulab.numpy.std(nums_list))
+timeit("ulab only, with ndarray", lambda: ulab.numpy.std(nums_array))
diff --git a/ulab_Crunch_Numbers_Fast/cluebarometer.py b/ulab_Crunch_Numbers_Fast/cluebarometer.py
index 3096b9a6b..9344a3461 100644
--- a/ulab_Crunch_Numbers_Fast/cluebarometer.py
+++ b/ulab_Crunch_Numbers_Fast/cluebarometer.py
@@ -3,8 +3,7 @@ import time
 import adafruit_bmp280
 import board
 import displayio
-import ulab
-import ulab.filter
+import ulab.numpy as np
 
 # Blank the screen.  Scrolling text causes unwanted delays.
 d = displayio.Group()
@@ -15,7 +14,7 @@ board.DISPLAY.show(d)
 # Transition bandwidth: 0.16Hz
 # Window type: Hamming
 # Filter has 311 coefficients
-taps = ulab.array([
+taps = np.array([
     -0.000050679794726066, -0.000041099278318167, -0.000031279920668665,
     -0.000021183486597150, -0.000010770285292045, +0.000000000000000000,
     +0.000011167446754809, +0.000022770999889941, +0.000034847558259864,
@@ -31,8 +30,8 @@ taps = ulab.array([
     +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.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,
@@ -63,8 +62,8 @@ taps = ulab.array([
     +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.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,
@@ -85,8 +84,8 @@ taps = ulab.array([
     +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.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,
@@ -143,7 +142,7 @@ 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))
+data = np.zeros(len(taps))
 t0 = deadline = time.monotonic_ns()
 n = 0
 # Take an initial reading to subtract off later, so that the graph in mu
@@ -163,9 +162,9 @@ while True:
         data = data + value
     else:
         # Otherwise, add it as the next sample
-        ulab.numerical.roll(data, 1)
+        data = np.roll(data, 1)
         data[-1] = value
-    filtered = ulab.numerical.sum(data * taps)
+    filtered = np.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
diff --git a/ulab_Crunch_Numbers_Fast/cluepulse.py b/ulab_Crunch_Numbers_Fast/cluepulse.py
index 36252a61e..3bf8beef4 100644
--- a/ulab_Crunch_Numbers_Fast/cluepulse.py
+++ b/ulab_Crunch_Numbers_Fast/cluepulse.py
@@ -3,8 +3,7 @@ import time
 import adafruit_apds9960.apds9960
 import board
 import digitalio
-import ulab
-import ulab.filter
+import ulab.numpy as np
 
 # Blank the screen.  Scrolling text causes unwanted delays.
 import displayio
@@ -18,7 +17,7 @@ board.DISPLAY.show(d)
 # Window type: Regular
 # Number of coefficients: 31
 # Manually trimmed to 16 coefficients
-taps = ulab.array([
+taps = np.array([
     +0.861745279666917052/2,
     -0.134728583242092248,
     -0.124472980501612152,
@@ -89,7 +88,7 @@ def main():
 
         # And our data structures
         # The most recent data samples, equal in number to the filter taps
-        data = ulab.zeros(len(taps))
+        data = np.zeros(len(taps))
         # The filtered value on the previous iteration
         old_value = 1
         # The times of the most recent pulses registered.  Increasing this number
@@ -111,11 +110,11 @@ def main():
             deadline += dt
             sleep_deadline(deadline)
             value = sum(sensor.color_data) # Combination of all channels
-            ulab.numerical.roll(data, 1)
+            data = np.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)
+            filtered = np.sum(data * taps)
 
             # We gathered enough data to fill the filters, and
             # the light value crossed the zero line in the positive direction
diff --git a/ulab_Crunch_Numbers_Fast/ledwave.py b/ulab_Crunch_Numbers_Fast/ledwave.py
index 54e838656..9dae735a5 100644
--- a/ulab_Crunch_Numbers_Fast/ledwave.py
+++ b/ulab_Crunch_Numbers_Fast/ledwave.py
@@ -2,9 +2,8 @@ import random
 
 import board
 import neopixel
-from _pixelbuf import wheel
-import ulab
-import ulab.filter
+from rainbowio import colorwheel as wheel
+import ulab.numpy as np
 
 # Customize your neopixel configuration here...
 pixel_pin = board.D5
@@ -12,11 +11,11 @@ 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.])
+ddt = np.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
+    deriv = np.convolve(u, ddt)[1:-1] * C2
     up = -um + u * 2 + deriv + f * dt2
     up[0] = 0
     up[n-1] = 0
@@ -29,11 +28,11 @@ def main():
     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)
+    u = np.zeros(num_pixels)
+    um = np.zeros(num_pixels)
+    f = np.zeros(num_pixels)
 
-    slope = ulab.linspace(0, 256, num=num_pixels)
+    slope = np.linspace(0, 256, num=num_pixels)
     th = 1
 
     # the first time is always random (is that a contradiction?)
diff --git a/ulab_Crunch_Numbers_Fast/waterfall.py b/ulab_Crunch_Numbers_Fast/waterfall.py
index 309839619..d34bebdb4 100644
--- a/ulab_Crunch_Numbers_Fast/waterfall.py
+++ b/ulab_Crunch_Numbers_Fast/waterfall.py
@@ -7,9 +7,8 @@ import array
 import board
 import audiobusio
 import displayio
-import ulab
-import ulab.fft
-import ulab.vector
+import ulab.numpy as np
+from ulab.scipy.signal import spectrogram
 
 display = board.DISPLAY
 
@@ -75,14 +74,14 @@ def main():
 
     while True:
         mic.record(samples_bit, len(samples_bit))
-        samples = ulab.array(samples_bit[3:])
-        spectrogram1 = ulab.fft.spectrogram(samples)
+        samples = np.array(samples_bit[3:])
+        spectrogram1 = spectrogram(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 = np.log(spectrogram1 + 1e-7)
         spectrogram1 = spectrogram1[1:(fft_size//2)-1]
-        min_curr = ulab.numerical.min(spectrogram1)
-        max_curr = ulab.numerical.max(spectrogram1)
+        min_curr = np.min(spectrogram1)
+        max_curr = np.max(spectrogram1)
 
         if max_curr > max_all:
             max_all = max_curr
@@ -94,7 +93,7 @@ def main():
         # Plot FFT
         data = (spectrogram1 - min_curr) * (51. / (max_all - min_curr))
         # This clamps any negative numbers to zero
-        data = data * ulab.array((data > 0))
+        data = data * np.array((data > 0))
         graph.show(data)
 
 main()