From fb1bd0b16917113274a5387e219175c330e344f8 Mon Sep 17 00:00:00 2001
From: Sasha Rush <srush.research@gmail.com>
Date: Fri, 16 Jul 2021 19:51:31 -0400
Subject: [PATCH 1/5] MFCC speech processing demo

---
 examples/alignment.dx | 381 ++++++++++++++++++++++++++++++++++++++++++
 1 file changed, 381 insertions(+)
 create mode 100644 examples/alignment.dx

diff --git a/examples/alignment.dx b/examples/alignment.dx
new file mode 100644
index 000000000..308f92937
--- /dev/null
+++ b/examples/alignment.dx
@@ -0,0 +1,381 @@
+import fft
+import parser
+import plot
+
+' # Speech Processing
+
+' This notebook describes the basic steps behind speech processing.
+ The goal will be to walk through the process from reading in a wave
+ file to being ready to train a full speech recognition system. 
+
+
+' ## Wave Files
+  To begin we need to be able to read in .wav files. To do this
+  we follow the wave file spec and define a header and the raw
+  data format.
+
+' [Wave Format Documentation](https://docs.fileformat.com/audio/wav/)
+  
+
+data WavHeader =
+     AsWavHeader {riff:(Fin 4=>Char) &
+                  size:Int &
+                  type:(Fin 4=>Char) &
+                  chunk:(Fin 4=>Char) &
+                  formatlength:Int &
+                  channels:Int & 
+                  samplerate:Int &
+                  bitspersample:Int &
+                  datasize:Int}
+                 
+data Wav = AsWav header:WavHeader dat:(List Int)
+
+' ### Binary Manipulation
+  We will need a couple additional binary conversion functions. (These are naive for now.)
+
+z = FToI (pow 2.0 15.0)
+def W32ToI (x : Word32): Int  =
+    y:Int = internalCast _ x
+    select (y <= z) y ((-1)*z + (y- z))
+       
+def W8ToW32 (x : Word8): Word32  = internalCast _  x
+
+def bytesToInt32 (chunk : Fin 4 => Byte) : Int =
+  [a, b, c, d] = map W8ToW32 chunk
+  W32ToI ((d .<<. 24) .|. (c .<<. 16) .|. (b .<<. 8) .|. a)
+
+
+def bytesToInt16 (chunk : Fin 2 => Byte) : Int =
+  [a, b] = map W8ToW32 chunk
+  W32ToI ((b .<<. 8) .|. a)
+
+
+' ### Parsing Types
+
+' We also add a couple additional parser helpers.
+
+def parseGroup (parser: Parser a) : (Parser (m => a)) = MkParser \h. for i. parse h parser
+parse4: Parser (Fin 4=>Char) = parseGroup parseAny
+parse2: Parser (Fin 2=>Char) = parseGroup parseAny
+
+def parse16 : Parser Int = MkParser \h.
+    bytesToInt16 $ parse h parse2
+
+def pString (x:String) : Parser Unit = MkParser \h.
+    (AsList n ls) = x
+    iter \i .
+         case i < n of
+              True ->
+                  _ = parse h (pChar ls.(i@_))
+                  Continue
+              False -> Done ()
+
+
+' ### Wave Parsing and IO
+
+' We define a parser to construct the header and to read in the data.
+
+wavparser : Parser Wav = MkParser \h.          
+              riff = parse h parse4
+              size = bytesToInt32 $ parse h parse4
+              type = parse h parse4
+              chunk = parse h parse4
+              _ = parse h parse16
+              formatlength = bytesToInt32 $ parse h parse4
+              channels = bytesToInt16 $ parse h parse2
+              samplerate = bytesToInt32 $ parse h parse4
+              _ = bytesToInt32 $ parse h parse4
+              _ = bytesToInt16 $ parse h parse2
+              bitspersample = bytesToInt16 $ parse h parse2
+              _ = parse h $ pString (AsList 4 ['d', 'a', 't', 'a'])
+              datasize = bytesToInt32 $ parse h parse4
+              x = parse h $ parseMany parse16
+              header = AsWavHeader {riff=riff, size=size,
+                     type=type, chunk=chunk,
+                     formatlength=formatlength,
+                     channels=channels, samplerate=samplerate,
+                     bitspersample=bitspersample,
+                     datasize=datasize}
+              AsWav header x
+
+
+' We add a helper function to read in the whole file. 
+
+def fullstream (stream : Stream ReadMode) : {IO} String =
+  yieldState (AsList _ []) \results.
+     iter \_. 
+           str = fread stream
+           (AsList n _) = str
+           results := (get results) <> str
+           case n == 0 of
+                True -> Done ()
+                False -> Continue
+                      
+x  = unsafeIO do runParserPartial (withFile "examples/speech2.wav" ReadMode fullstream) wavparser
+(AsList len raw_sample_wav) = case x of
+     Nothing -> mempty
+     (Just (AsWav header dat)) -> dat
+
+Datsize = Fin len
+samplerate = 8000.0
+
+
+' ### The Audio
+
+
+' We are going to work with an MNist style of audio that has short snippets of
+  people saying letters.
+
+s = unsafeIO do base64Encode (withFile "examples/speech2.wav" ReadMode fullstream)
+
+:html "<audio controls src='data:audio/wav;base64,"<> s <>"'/>"
+
+
+' Here is what this looks like as a wave. 
+
+:html showPlot $ xyPlot (for i. IToF (ordinal i))  (for i. IToF raw_sample_wav.i)
+
+
+' ## MFCC Signal Processing
+
+' Before doing machine learning on the data we need to process it.
+  We are going to follow a standard pipeline to produce MFCC features
+  for the data. The steps for this process are outlined in wikipedia.
+  
+'  https://en.wikipedia.org/wiki/Mel-frequency_cepstrum
+
+' ### Signal Tools
+
+' First we will define an unsafe window function that will let us easily take convolutions.
+
+def window [Add out] (stride:Int) (ls: m => out) : (n => window => out) =
+    window_size = size window
+    for i: n.
+        case ((ordinal i) * stride + window_size) < (size m) of
+            True ->
+               j = ((ordinal i) * stride) @ m
+               k = ((ordinal i) * stride + window_size) @ m
+               x = for l':(j..k).
+                  l = %inject l'
+                  ls.l
+               unsafeCastTable window x
+ 
+            False -> zero
+
+' Add arbitrary padding based on type.
+
+def pad [Add out] (tab: m => out) : (n => out) =
+    for i.
+        case ordinal i < size m of
+             True -> tab.((ordinal i)@_)
+             False -> zero
+
+
+' ### Step 0: Preemphasis Smoothing
+
+' As a warmup, we will play around with windowing by running a simple
+convolution filter over the data. This preemphasis filter smooths out the signal.
+
+-- Hyperparam
+smoothing_coef = 0.97
+def smoothing (x:Fin 2=>Float)  : Float = x.(1@_) - (smoothing_coef * x.(0@_))
+    
+signal : Datsize => Float = map smoothing $ window 1 (for i. IToF raw_sample_wav.i)
+
+
+' ### Step 1: Movement to Frequency domain
+
+' Next we move from time to frequency domain by applying an FFT.
+
+' Instead of applying on FFT on the whole sequence we break it up into overlapping windows and apply an FFT on each one.  The windows in speech are known as frames, they are typically 0.025 seconds long. We use a stride of 0.01 seconds for overlap.
+
+-- Hyperparam
+WindowSize = Fin (FToI (0.025 * samplerate))
+stride = (FToI (0.01 * samplerate))
+
+
+PostWindow = Fin (idiv (size Datsize) stride)
+
+frame_split : PostWindow => WindowSize => Float = window stride signal
+
+
+' One issue though is that the FFT does not change the size of its input,
+  so if we have a short sequence we will end up with low frequency space
+  resolution. We get around this by padding out the time sequence.
+
+NFFT = Fin 512
+
+' Finally we only keep the positive components of the FFT which are the first half
+  plus one.
+
+Positive = Fin ((idiv 512  2) + 1)
+
+fft_frames = for i.
+    prefft: NFFT => Float = pad frame_split.i
+    ff = fft_real prefft
+    for j:Positive. (complex_abs ff.((ordinal j)@NFFT)) / (IToF (size NFFT))
+
+
+:html matshow fft_frames
+
+
+' ## Step 2 / 3: Construct Scales
+
+' The next problem is that we have the data in frequency-space which is not that useful for the
+  problem of speech recognition. We would like to move the data to a scale that is more focused on
+  the range of human communication and our biology.
+
+
+' ### Mel FilterBank
+
+' To do this we use a nonlinear rescaling known as the Mel Scale. 
+
+' Mel scale looks like this.
+
+Hertz = Float
+Mel = Float
+def mel (f:Hertz) : Mel = 1125.0 * log (1.0 + (f / 700.0))
+def imel (f:Mel) : Hertz = 700.0 * (exp (f / 1125.0) - 1.0)
+
+' However it is now so easy to move between frequency scales. To do this we will need
+' a bit of Dex machinery.
+
+' ### Linear Scales
+
+' Let us define a scaled range as a linear mapping between ordinals and a range. 
+
+LinearScale = {start:Float & end:Float}
+data ScaledRange a = AsScaledRange scale:LinearScale
+
+' Given a scale, we can easily map from an index to a point, get the step, get the full domain,
+  and take the binned inverse to get the nearest index.
+
+def step (AsScaledRange {start, end}:ScaledRange a) : Float =
+    (1.0 / (IToF ((size a) - 1) )) * (end - start)
+def scaled (AsScaledRange {start, end}:ScaledRange a) (x:a): Float =
+    start + ((IToF (ordinal x)) / (IToF ((size a) - 1))) * (end - start)
+def domain (scale:ScaledRange a) : a => Float =
+    for i. scaled scale i
+def bin (scale:ScaledRange a) (x:Float) : a =
+    (AsScaledRange {start, end}) = scale
+    (FToI (floor ((x - start) / (step scale))))@_
+
+' ## Mapping Between Scales
+
+' We can then define a two linear ranges. One in Hertz (corresponding to the FFT) and one in Mel Space.
+
+MelBins = Fin 28
+
+hscale : ScaledRange Positive = AsScaledRange {start=0.0, end=samplerate / 2.0}
+melscale : ScaledRange MelBins = AsScaledRange {start=mel 0.0, end=mel (samplerate / 2.0)}
+domain melscale
+
+mel_to_hscale = for i: MelBins.
+        mel = scaled melscale i
+        bin hscale (imel mel)
+
+
+' ## Triangle Filter Bank
+
+' In order to apply this mapping we need a way to handle the non-linear disconnect. To do this
+  We build a triangle windowing function.
+
+:t fft_frames
+
+def triangle_window (mpt:m) (mlower:m) (mupper:m) (scale: ScaledRange m) : m => Float =
+    pt = scaled scale mpt
+    lower = scaled scale mlower
+    upper = scaled scale mupper
+    for i.
+        k = scaled scale i
+        if (k >= lower) && (k <= pt) 
+           then (k - lower) / (pt - lower)
+           else if (k <= upper) && (k >= pt) 
+             then (upper - k) / (upper - pt)
+             else 0.0
+
+' The 25'th Mel bin is going to take in the mass from this range of hertz's bins.
+
+:html showPlot $ xyPlot (domain hscale) (triangle_window mel_to_hscale.(25@_) mel_to_hscale.(24@_) mel_to_hscale.(26@_) hscale)
+
+' Whereas the 19'th Mel bin is going to take in the mass from this range of hertz's bins.
+
+:html showPlot $ xyPlot (domain hscale) (triangle_window mel_to_hscale.(19@_) mel_to_hscale.(18@_) mel_to_hscale.(20@_) hscale)
+
+
+' Here is the full matrix.
+
+start = 1@MelBins
+end = 27@MelBins
+MelBinsIn = Fin 26
+H : MelBinsIn => Positive => Float = unsafeCastTable MelBinsIn $ for m' : (start..<end).
+      m = %inject m'
+      m1 = ((ordinal m) - 1)@MelBins
+      m2 = ((ordinal m) + 1)@MelBins
+      triangle_window mel_to_hscale.m mel_to_hscale.m1 mel_to_hscale.m2 hscale
+
+
+' We can plot these over each other to get the triangles for each Mel bin.
+
+d = for i. plotToDiagram $ xyPlot (domain hscale) H.i
+:html renderScaledSVG $ concatDiagrams d
+
+' The final Mel features are constructed applying a matrix multiply.
+
+eps = 0.000001
+melled : PostWindow => MelBinsIn => Float = for i. for m.  log ((sum for j. fft_frames.i.j * H.m.j) + eps)
+:html matshow for i j. melled.j.i
+
+' ## Step 4: Transform the Mel spectrum
+
+' The last step is to take a Discrete cosine transform in feature space. This will allow us to
+  further compress the signal down for usage. We start with 26 Mel features and after this step
+  we will end with 13 MFCCs. 
+
+' There is a nice property that the discrete cosine transform can be computed by running a
+  standard FFT on a transformed version of the data. We implement this based on 
+
+' [Fast Cosine Transform](https://dsp.stackexchange.com/questions/2807/fast-cosine-transform-via-fft)
+
+AB = {A : Unit | B : Unit}
+     
+def dct (input : m=>Float) : m=>Complex =
+   d = for (k, i, j) : (AB & m & AB). 
+        case j of
+             {|B=()|} -> case k of
+               {|A=()|} -> input.i
+               {|B=()|} -> input.(((size m)-1-(ordinal i))@_)
+             {|A=()|} -> 0.0
+   d' = unsafeCastTable (AB & AB & m) $ fft_real d
+   for i. d'.({|A=()|}, {|A=()|}, i)
+
+' Now we just apply DCT along each Mel window to get our final values.
+From this we drop the top DCT coefficient and keep the next 13.
+
+MFCCStart = 1@MelBinsIn  
+MFCCLen = 13@MelBinsIn  
+mfcc = for i.
+    d = dct melled.i
+    for j': (MFCCStart..MFCCLen).
+        j = %inject j'
+        (MkComplex m _) =  d.j
+        m
+
+' Finally there is one last heuristic (lifting) that is applied. This readjusts the
+  coefficients based on higher-frequencies. 
+
+lift_coef = 22.0
+
+mfcc_lifted = for i j.
+        -- Lifting
+        n = IToF (ordinal j)
+        lift = 1.0 + (lift_coef/2.)*sin(pi* n /lift_coef)
+        mfcc.i.j * lift
+
+
+:html matshow $ for i j. mfcc_lifted.j.i
+
+
+' We now have our MFCC features which are ready for training or analysing speech data.
+
+

From b77105f89274b4ac38509decf5bf780e93354015 Mon Sep 17 00:00:00 2001
From: Sasha Rush <srush.research@gmail.com>
Date: Fri, 16 Jul 2021 19:52:44 -0400
Subject: [PATCH 2/5] fft

---
 lib/fft.dx | 150 +++++++++++++++++++++++++++++++++++++++++++++++++++++
 1 file changed, 150 insertions(+)
 create mode 100644 lib/fft.dx

diff --git a/lib/fft.dx b/lib/fft.dx
new file mode 100644
index 000000000..5dd68c231
--- /dev/null
+++ b/lib/fft.dx
@@ -0,0 +1,150 @@
+'# Fast Fourier Transform
+For arrays whose size is a power of 2, we use a radix-2 algorithm based
+on the [Fhutark demo](https://github.com/diku-dk/fft/blob/master/lib/github.com/diku-dk/fft/stockham-radix-2.fut#L30).
+For non-power-of-2 sized arrays, it uses
+[Bluestein's Algorithm](https://en.wikipedia.org/wiki/Chirp_Z-transform),
+which calls the power-of-2 FFT as a subroutine.
+
+
+'### Helper functions
+
+def odd_sized_palindrome (mid:a) (seq:n=>a) :
+  ({backward:n | mid:Unit | zforward:n}=>a) =  -- Alphabetical order matters here.
+  -- Turns sequence 12345 into 543212345.
+  for i.
+    case i of
+      {|backward=i|} -> seq.(reflect i)
+      {|mid=i|} -> mid
+      {|zforward=i|} -> seq.i
+
+
+'## Inner FFT functions
+
+data FTDirection =
+  ForwardFT
+  InverseFT
+
+def butterfly_ixs (j':halfn) (pow2:Int) : (n & n & n & n) =
+  -- Re-index at a finer frequency.
+  -- halfn must have half the size of n.
+  -- For explanation, see https://en.wikipedia.org/wiki/Butterfly_diagram
+  -- Note: with fancier index sets, this might be replacable by reshapes.
+  j = ordinal j'
+  k = ((idiv j pow2) * pow2 * 2) + mod j pow2
+  left_write_ix  = unsafeFromOrdinal n k
+  right_write_ix = unsafeFromOrdinal n (k + pow2)
+
+  left_read_ix  = unsafeFromOrdinal n j
+  right_read_ix = unsafeFromOrdinal n (j + size halfn)
+  (left_read_ix, right_read_ix, left_write_ix, right_write_ix)
+
+def power_of_2_fft (direction: FTDirection) (x: n=>Complex) : n=>Complex =
+  -- Input size must be a power of 2.
+  -- Can enforce this with tables-as-index-sets like:
+  -- (x: (log2n=>(Fin 2))=>Complex)) once that's supported.
+  dir_const = case direction of
+    ForwardFT -> -pi
+    InverseFT -> pi
+
+  log2n = intlog2 (size n)
+  halfn = idiv (size n) 2
+  
+  ans = yieldState x \xRef.
+    for i:(Fin log2n).
+      ipow2 = intpow 2 (ordinal i)
+      xRef := yieldAccum (AddMonoid Complex) \bufRef.
+        for j:(Fin halfn).  -- Executes in parallel.
+          (left_read_ix, right_read_ix,
+           left_write_ix, right_write_ix) = butterfly_ixs j ipow2
+
+          -- Read one element from the last buffer, scaled.
+          angle = dir_const * (IToF $ mod (ordinal j) ipow2) / IToF ipow2
+          v = (get xRef!right_read_ix) * (MkComplex (cos angle) (sin angle))
+
+          -- Add and subtract it to the relevant places in the new buffer.
+          bufRef!left_write_ix  += (get (xRef!left_read_ix)) + v
+          bufRef!right_write_ix += (get (xRef!left_read_ix)) - v
+
+  case direction of
+    ForwardFT -> ans
+    InverseFT -> ans / (IToF (size n))
+
+def convolve_complex (u:n=>Complex) (v:m=>Complex) :
+  ({orig_vals:n | padding:m }=>Complex) =  -- Alphabetical order matters here.
+  -- Convolve by pointwise multiplication in the Fourier domain.
+  convolved_size = (size n) + (size m) - 1
+  working_size = nextpow2 convolved_size
+  u_padded = padTo (Fin working_size) zero u
+  v_padded = padTo (Fin working_size) zero v
+  spectral_u = power_of_2_fft ForwardFT u_padded
+  spectral_v = power_of_2_fft ForwardFT v_padded
+  spectral_conv = for i. spectral_u.i * spectral_v.i
+  padded_conv = power_of_2_fft InverseFT spectral_conv
+  slice padded_conv 0 {orig_vals:n | padding:m }
+
+def convolve (u:n=>Float) (v:m=>Float) : ({orig_vals:n | padding:m }=>Float) =
+  u' = for i. MkComplex u.i 0.0
+  v' = for i. MkComplex v.i 0.0
+  ans = convolve_complex u' v'
+  for i.
+    (MkComplex real imag) = ans.i
+    real
+
+
+'## FFT Interface
+
+def fft (x: n=>Complex): n=>Complex =
+  if isPowerOf2 (size n)
+    then power_of_2_fft ForwardFT x
+    else
+      -- Bluestein's algorithm.
+      -- Converts the general FFT into a convolution,
+      -- which is then solved with calls to a power-of-2 FFT.
+      im = MkComplex 0.0 1.0
+      wks = for i.
+        i_squared = IToF $ sq $ ordinal i
+        exp $ (-im) * (MkComplex (pi * i_squared / (IToF (size n))) 0.0)
+
+      (AsList _ tailTable) = tail wks 1
+      back_and_forth = odd_sized_palindrome (head wks) tailTable
+      xq = for i. x.i * wks.i
+      back_and_forth_conj = for i. complex_conj back_and_forth.i
+      convolution = convolve_complex xq back_and_forth_conj
+      convslice = slice convolution (size n - 1) n
+      for i. wks.i * convslice.i
+
+def ifft (xs: n=>Complex): n=>Complex =
+  if isPowerOf2 (size n)
+    then power_of_2_fft InverseFT xs
+    else
+      unscaled_fft = fft (for i. complex_conj xs.i)
+      for i. (complex_conj unscaled_fft.i) / (IToF (size n))
+
+def  fft_real (x: n=>Float): n=>Complex =  fft for i. MkComplex x.i 0.0
+def ifft_real (x: n=>Float): n=>Complex = ifft for i. MkComplex x.i 0.0
+
+def fft2 (x: n=>m=>Complex): n=>m=>Complex =
+  x'      = for i. fft x.i
+  transpose for i. fft (transpose x').i
+
+def ifft2 (x: n=>m=>Complex): n=>m=>Complex =
+  x'      = for i. ifft x.i
+  transpose for i. ifft (transpose x').i
+
+def  fft2_real (x: n=>m=>Float): n=>m=>Complex =  fft2 for i j. MkComplex x.i.j 0.0
+def ifft2_real (x: n=>m=>Float): n=>m=>Complex = ifft2 for i j. MkComplex x.i.j 0.0
+
+-------- Tests --------
+
+a = for i. MkComplex [10.1, -2.2, 8.3, 4.5, 9.3].i 0.0
+b = for i:(Fin 20) j:(Fin 70).
+  MkComplex (randn $ ixkey2 (newKey 0) i j) 0.0
+
+:p a ~~ (ifft $ fft a)
+> True
+:p a ~~ (fft $ ifft a)
+> True
+:p b ~~ (ifft2 $ fft2 b)
+> True
+:p b ~~ (fft2 $ ifft2 b)
+> True

From 26c9d8b6fdd77b29db38b89a506d458241892936 Mon Sep 17 00:00:00 2001
From: Sasha Rush <srush.research@gmail.com>
Date: Fri, 16 Jul 2021 20:01:36 -0400
Subject: [PATCH 3/5] scael

---
 examples/alignment.dx | 18 +++++++++---------
 1 file changed, 9 insertions(+), 9 deletions(-)

diff --git a/examples/alignment.dx b/examples/alignment.dx
index 308f92937..d9b5c7901 100644
--- a/examples/alignment.dx
+++ b/examples/alignment.dx
@@ -283,15 +283,15 @@ mel_to_hscale = for i: MelBins.
 :t fft_frames
 
 def triangle_window (mpt:m) (mlower:m) (mupper:m) (scale: ScaledRange m) : m => Float =
-    pt = scaled scale mpt
-    lower = scaled scale mlower
-    upper = scaled scale mupper
-    for i.
-        k = scaled scale i
-        if (k >= lower) && (k <= pt) 
-           then (k - lower) / (pt - lower)
-           else if (k <= upper) && (k >= pt) 
-             then (upper - k) / (upper - pt)
+    pt = IToF (ordinal mpt)
+    lower = IToF (ordinal mlower)
+    upper = IToF (ordinal mupper)
+    for i'.
+        i = IToF (ordinal i')
+        if (i >= lower) && (i <= pt) 
+           then (i - lower) / (pt - lower)
+           else if (i <= upper) && (i >= pt) 
+             then (upper - i) / (upper - pt)
              else 0.0
 
 ' The 25'th Mel bin is going to take in the mass from this range of hertz's bins.

From 31317d6d11f3936b25d9e029e9e644fb1764193a Mon Sep 17 00:00:00 2001
From: Sasha Rush <srush.research@gmail.com>
Date: Wed, 21 Jul 2021 00:25:57 -0400
Subject: [PATCH 4/5] Move around some of the functions

---
 examples/alignment.dx | 228 +++++++++++++++++-------------------------
 examples/fft.dx       | 150 ---------------------------
 lib/parser.dx         |   6 ++
 lib/prelude.dx        |  62 +++++++++++-
 4 files changed, 161 insertions(+), 285 deletions(-)
 delete mode 100644 examples/fft.dx

diff --git a/examples/alignment.dx b/examples/alignment.dx
index d9b5c7901..c4257d548 100644
--- a/examples/alignment.dx
+++ b/examples/alignment.dx
@@ -6,7 +6,7 @@ import plot
 
 ' This notebook describes the basic steps behind speech processing.
  The goal will be to walk through the process from reading in a wave
- file to being ready to train a full speech recognition system. 
+ file to being ready to train a full speech recognition system.
 
 
 ' ## Wave Files
@@ -15,19 +15,18 @@ import plot
   data format.
 
 ' [Wave Format Documentation](https://docs.fileformat.com/audio/wav/)
-  
+
 
 data WavHeader =
-     AsWavHeader {riff:(Fin 4=>Char) &
-                  size:Int &
+     AsWavHeader {size:Int &
                   type:(Fin 4=>Char) &
                   chunk:(Fin 4=>Char) &
                   formatlength:Int &
-                  channels:Int & 
+                  channels:Int &
                   samplerate:Int &
                   bitspersample:Int &
                   datasize:Int}
-                 
+
 data Wav = AsWav header:WavHeader dat:(List Int)
 
 ' ### Binary Manipulation
@@ -37,7 +36,7 @@ z = FToI (pow 2.0 15.0)
 def W32ToI (x : Word32): Int  =
     y:Int = internalCast _ x
     select (y <= z) y ((-1)*z + (y- z))
-       
+
 def W8ToW32 (x : Word8): Word32  = internalCast _  x
 
 def bytesToInt32 (chunk : Fin 4 => Byte) : Int =
@@ -54,84 +53,59 @@ def bytesToInt16 (chunk : Fin 2 => Byte) : Int =
 
 ' We also add a couple additional parser helpers.
 
-def parseGroup (parser: Parser a) : (Parser (m => a)) = MkParser \h. for i. parse h parser
-parse4: Parser (Fin 4=>Char) = parseGroup parseAny
-parse2: Parser (Fin 2=>Char) = parseGroup parseAny
-
-def parse16 : Parser Int = MkParser \h.
-    bytesToInt16 $ parse h parse2
-
-def pString (x:String) : Parser Unit = MkParser \h.
-    (AsList n ls) = x
-    iter \i .
-         case i < n of
-              True ->
-                  _ = parse h (pChar ls.(i@_))
-                  Continue
-              False -> Done ()
+def parseChars (n:Int) : (Parser (Fin n => Char)) =
+    MkParser \h. for i. parse h parseAny
 
+def parseI16 : Parser Int = MkParser \h.
+    bytesToInt16 $ parse h $ parseChars 2
 
 ' ### Wave Parsing and IO
 
 ' We define a parser to construct the header and to read in the data.
 
-wavparser : Parser Wav = MkParser \h.          
-              riff = parse h parse4
-              size = bytesToInt32 $ parse h parse4
-              type = parse h parse4
-              chunk = parse h parse4
-              _ = parse h parse16
-              formatlength = bytesToInt32 $ parse h parse4
-              channels = bytesToInt16 $ parse h parse2
-              samplerate = bytesToInt32 $ parse h parse4
-              _ = bytesToInt32 $ parse h parse4
-              _ = bytesToInt16 $ parse h parse2
-              bitspersample = bytesToInt16 $ parse h parse2
-              _ = parse h $ pString (AsList 4 ['d', 'a', 't', 'a'])
-              datasize = bytesToInt32 $ parse h parse4
-              x = parse h $ parseMany parse16
-              header = AsWavHeader {riff=riff, size=size,
-                     type=type, chunk=chunk,
-                     formatlength=formatlength,
-                     channels=channels, samplerate=samplerate,
-                     bitspersample=bitspersample,
-                     datasize=datasize}
+wavparser : Parser Wav = MkParser \h.
+              parse h $ parseChars 4
+              size = bytesToInt32 $ parse h $ parseChars 4
+              type = parse h $ parseChars 4
+              chunk = parse h $ parseChars 4
+              parse h parseI16
+              formatlength = bytesToInt32 $ parse h $ parseChars 4
+              channels = parse h parseI16
+              samplerate = bytesToInt32 $ parse h $ parseChars 4
+              bytesToInt32 $ parse h $ parseChars 4
+              parse h parseI16
+              bitspersample = parse h parseI16
+              parse h $ pString (AsList 4 ['d', 'a', 't', 'a'])
+              datasize = bytesToInt32 $ parse h $ parseChars 4
+              x = parse h $ parseMany parseI16
+              header = AsWavHeader {size,
+                     type, chunk,
+                     formatlength,
+                     channels, samplerate,
+                     bitspersample,
+                     datasize}
               AsWav header x
 
 
-' We add a helper function to read in the whole file. 
-
-def fullstream (stream : Stream ReadMode) : {IO} String =
-  yieldState (AsList _ []) \results.
-     iter \_. 
-           str = fread stream
-           (AsList n _) = str
-           results := (get results) <> str
-           case n == 0 of
-                True -> Done ()
-                False -> Continue
-                      
-x  = unsafeIO do runParserPartial (withFile "examples/speech2.wav" ReadMode fullstream) wavparser
+x  = unsafeIO do runParserPartial (readFile "examples/speech2.wav") wavparser
 (AsList len raw_sample_wav) = case x of
      Nothing -> mempty
      (Just (AsWav header dat)) -> dat
 
-Datsize = Fin len
-samplerate = 8000.0
 
+samplerate = 8000
 
 ' ### The Audio
 
-
 ' We are going to work with an MNist style of audio that has short snippets of
   people saying letters.
 
-s = unsafeIO do base64Encode (withFile "examples/speech2.wav" ReadMode fullstream)
+s = unsafeIO do base64Encode (readFile "examples/speech2.wav")
 
 :html "<audio controls src='data:audio/wav;base64,"<> s <>"'/>"
 
 
-' Here is what this looks like as a wave. 
+' Here is what this looks like as a wave.
 
 :html showPlot $ xyPlot (for i. IToF (ordinal i))  (for i. IToF raw_sample_wav.i)
 
@@ -141,47 +115,31 @@ s = unsafeIO do base64Encode (withFile "examples/speech2.wav" ReadMode fullstrea
 ' Before doing machine learning on the data we need to process it.
   We are going to follow a standard pipeline to produce MFCC features
   for the data. The steps for this process are outlined in wikipedia.
-  
-'  https://en.wikipedia.org/wiki/Mel-frequency_cepstrum
 
-' ### Signal Tools
+' [MFCC Wikipedia](https://en.wikipedia.org/wiki/Mel-frequency_cepstrum)
 
-' First we will define an unsafe window function that will let us easily take convolutions.
+' We will rely pretty heavily on the windowing tools provided in the prelude.
 
-def window [Add out] (stride:Int) (ls: m => out) : (n => window => out) =
-    window_size = size window
-    for i: n.
-        case ((ordinal i) * stride + window_size) < (size m) of
-            True ->
-               j = ((ordinal i) * stride) @ m
-               k = ((ordinal i) * stride + window_size) @ m
-               x = for l':(j..k).
-                  l = %inject l'
-                  ls.l
-               unsafeCastTable window x
- 
-            False -> zero
+' ### Step 0: Preemphasis Smoothing
 
-' Add arbitrary padding based on type.
 
-def pad [Add out] (tab: m => out) : (n => out) =
-    for i.
-        case ordinal i < size m of
-             True -> tab.((ordinal i)@_)
-             False -> zero
+' We begin by padding so it is divisible by 80. This will allow us
+  to keep the sizes throughout.
 
-
-' ### Step 0: Preemphasis Smoothing
+Datsize = (Pad (Fin 0) (Fin len) (Fin (80 -(mod len 80))))
+sample_wav : Datsize => Float =
+           pad 0.0 (for j. IToF raw_sample_wav.j)
 
 ' As a warmup, we will play around with windowing by running a simple
 convolution filter over the data. This preemphasis filter smooths out the signal.
 
+
 -- Hyperparam
 smoothing_coef = 0.97
-def smoothing (x:Fin 2=>Float)  : Float = x.(1@_) - (smoothing_coef * x.(0@_))
-    
-signal : Datsize => Float = map smoothing $ window 1 (for i. IToF raw_sample_wav.i)
+def smoothing (x:(Window Unit (Fin 0))=>Float)  : Float =
+    x.{|mid=()|} - (smoothing_coef * x.{|left=()|})
 
+signal : Datsize => Float = map smoothing $ window $ pad 0.0 sample_wav
 
 ' ### Step 1: Movement to Frequency domain
 
@@ -190,13 +148,12 @@ signal : Datsize => Float = map smoothing $ window 1 (for i. IToF raw_sample_wav
 ' Instead of applying on FFT on the whole sequence we break it up into overlapping windows and apply an FFT on each one.  The windows in speech are known as frames, they are typically 0.025 seconds long. We use a stride of 0.01 seconds for overlap.
 
 -- Hyperparam
-WindowSize = Fin (FToI (0.025 * samplerate))
-stride = (FToI (0.01 * samplerate))
+FrameWindow = Window (Fin 0) (Fin 199)
+Step = Fin  $ idiv samplerate 80
+PostWindow = Fin $ idiv (size Datsize) (size Step)
 
-
-PostWindow = Fin (idiv (size Datsize) stride)
-
-frame_split : PostWindow => WindowSize => Float = window stride signal
+frame_split : PostWindow => FrameWindow => Float =
+             stride $ castTable (_ & Step) $ window $ pad 0.0 signal
 
 
 ' One issue though is that the FFT does not change the size of its input,
@@ -209,11 +166,12 @@ NFFT = Fin 512
   plus one.
 
 Positive = Fin ((idiv 512  2) + 1)
+PaddingFFT = Fin ((size NFFT) - (size FrameWindow))
 
 fft_frames = for i.
-    prefft: NFFT => Float = pad frame_split.i
+    prefft : (Pad (Fin 0) FrameWindow PaddingFFT) => Float  = pad 0.0 frame_split.i
     ff = fft_real prefft
-    for j:Positive. (complex_abs ff.((ordinal j)@NFFT)) / (IToF (size NFFT))
+    for j:Positive. (complex_abs ff.((ordinal j)@_)) / (IToF (size NFFT))
 
 
 :html matshow fft_frames
@@ -228,7 +186,7 @@ fft_frames = for i.
 
 ' ### Mel FilterBank
 
-' To do this we use a nonlinear rescaling known as the Mel Scale. 
+' To do this we use a nonlinear rescaling known as the Mel Scale.
 
 ' Mel scale looks like this.
 
@@ -242,7 +200,7 @@ def imel (f:Mel) : Hertz = 700.0 * (exp (f / 1125.0) - 1.0)
 
 ' ### Linear Scales
 
-' Let us define a scaled range as a linear mapping between ordinals and a range. 
+' Let us define a scaled range as a linear mapping between ordinals and a range. We keep the type around so we know the size of the domain.
 
 LinearScale = {start:Float & end:Float}
 data ScaledRange a = AsScaledRange scale:LinearScale
@@ -264,13 +222,15 @@ def bin (scale:ScaledRange a) (x:Float) : a =
 
 ' We can then define a two linear ranges. One in Hertz (corresponding to the FFT) and one in Mel Space.
 
-MelBins = Fin 28
+MelBins = Fin 26
+
 
-hscale : ScaledRange Positive = AsScaledRange {start=0.0, end=samplerate / 2.0}
-melscale : ScaledRange MelBins = AsScaledRange {start=mel 0.0, end=mel (samplerate / 2.0)}
+top = (IToF samplerate) / 2.0
+hscale : ScaledRange Positive = AsScaledRange {start=0.0, end=top}
+melscale : ScaledRange (Pad Unit MelBins Unit) = AsScaledRange {start=mel 0.0, end=mel top}
 domain melscale
 
-mel_to_hscale = for i: MelBins.
+mel_to_hscale = for i.
         mel = scaled melscale i
         bin hscale (imel mel)
 
@@ -282,37 +242,38 @@ mel_to_hscale = for i: MelBins.
 
 :t fft_frames
 
-def triangle_window (mpt:m) (mlower:m) (mupper:m) (scale: ScaledRange m) : m => Float =
-    pt = IToF (ordinal mpt)
-    lower = IToF (ordinal mlower)
-    upper = IToF (ordinal mupper)
-    for i'.
-        i = IToF (ordinal i')
-        if (i >= lower) && (i <= pt) 
-           then (i - lower) / (pt - lower)
-           else if (i <= upper) && (i >= pt) 
-             then (upper - i) / (upper - pt)
-             else 0.0
+def triangle_window (mlower:m) (mpt:m)  (mupper:m) : m => Float =
+    up = (mlower..mpt)
+    down = (mpt..mupper)
+    uprange = AsScaledRange {start=0.0, end=1.0}
+    downrange = AsScaledRange {start=1.0, end=0.0}
+    yieldState zero \ref.
+        for i' : up.
+            i = %inject i'
+            (ref!i) := scaled uprange i'
+        for i' : down.
+            i = %inject i'
+            (ref!i) := scaled downrange i'
+
 
 ' The 25'th Mel bin is going to take in the mass from this range of hertz's bins.
 
-:html showPlot $ xyPlot (domain hscale) (triangle_window mel_to_hscale.(25@_) mel_to_hscale.(24@_) mel_to_hscale.(26@_) hscale)
+:html showPlot $ xyPlot (domain hscale) (triangle_window mel_to_hscale.(24@_) mel_to_hscale.(25@_) mel_to_hscale.(26@_))
 
 ' Whereas the 19'th Mel bin is going to take in the mass from this range of hertz's bins.
 
-:html showPlot $ xyPlot (domain hscale) (triangle_window mel_to_hscale.(19@_) mel_to_hscale.(18@_) mel_to_hscale.(20@_) hscale)
+:html showPlot $ xyPlot (domain hscale) (triangle_window mel_to_hscale.(18@_) mel_to_hscale.(19@_) mel_to_hscale.(20@_))
 
 
 ' Here is the full matrix.
 
-start = 1@MelBins
-end = 27@MelBins
-MelBinsIn = Fin 26
-H : MelBinsIn => Positive => Float = unsafeCastTable MelBinsIn $ for m' : (start..<end).
-      m = %inject m'
-      m1 = ((ordinal m) - 1)@MelBins
-      m2 = ((ordinal m) + 1)@MelBins
-      triangle_window mel_to_hscale.m mel_to_hscale.m1 mel_to_hscale.m2 hscale
+
+H : MelBins => Positive => Float =
+  map (\x. triangle_window x.{|left=()|}  x.{|mid=()|} x.{|right=()|}) $
+     window $ castTable (Pad Unit MelBins Unit) mel_to_hscale
+
+
+
 
 
 ' We can plot these over each other to get the triangles for each Mel bin.
@@ -323,24 +284,25 @@ d = for i. plotToDiagram $ xyPlot (domain hscale) H.i
 ' The final Mel features are constructed applying a matrix multiply.
 
 eps = 0.000001
-melled : PostWindow => MelBinsIn => Float = for i. for m.  log ((sum for j. fft_frames.i.j * H.m.j) + eps)
+melled : PostWindow => MelBins => Float =
+       for i. for m. log ((sum for j. fft_frames.i.j * H.m.j) + eps)
 :html matshow for i j. melled.j.i
 
 ' ## Step 4: Transform the Mel spectrum
 
 ' The last step is to take a Discrete cosine transform in feature space. This will allow us to
   further compress the signal down for usage. We start with 26 Mel features and after this step
-  we will end with 13 MFCCs. 
+  we will end with 13 MFCCs.
 
 ' There is a nice property that the discrete cosine transform can be computed by running a
-  standard FFT on a transformed version of the data. We implement this based on 
+  standard FFT on a transformed version of the data. We implement this based on
 
 ' [Fast Cosine Transform](https://dsp.stackexchange.com/questions/2807/fast-cosine-transform-via-fft)
 
 AB = {A : Unit | B : Unit}
-     
+
 def dct (input : m=>Float) : m=>Complex =
-   d = for (k, i, j) : (AB & m & AB). 
+   d = for (k, i, j) : (AB & m & AB).
         case j of
              {|B=()|} -> case k of
                {|A=()|} -> input.i
@@ -352,8 +314,8 @@ def dct (input : m=>Float) : m=>Complex =
 ' Now we just apply DCT along each Mel window to get our final values.
 From this we drop the top DCT coefficient and keep the next 13.
 
-MFCCStart = 1@MelBinsIn  
-MFCCLen = 13@MelBinsIn  
+MFCCStart = 1@MelBins
+MFCCLen = 13@MelBins
 mfcc = for i.
     d = dct melled.i
     for j': (MFCCStart..MFCCLen).
@@ -362,7 +324,7 @@ mfcc = for i.
         m
 
 ' Finally there is one last heuristic (lifting) that is applied. This readjusts the
-  coefficients based on higher-frequencies. 
+  coefficients based on higher-frequencies.
 
 lift_coef = 22.0
 
@@ -377,5 +339,3 @@ mfcc_lifted = for i j.
 
 
 ' We now have our MFCC features which are ready for training or analysing speech data.
-
-
diff --git a/examples/fft.dx b/examples/fft.dx
deleted file mode 100644
index 5dd68c231..000000000
--- a/examples/fft.dx
+++ /dev/null
@@ -1,150 +0,0 @@
-'# Fast Fourier Transform
-For arrays whose size is a power of 2, we use a radix-2 algorithm based
-on the [Fhutark demo](https://github.com/diku-dk/fft/blob/master/lib/github.com/diku-dk/fft/stockham-radix-2.fut#L30).
-For non-power-of-2 sized arrays, it uses
-[Bluestein's Algorithm](https://en.wikipedia.org/wiki/Chirp_Z-transform),
-which calls the power-of-2 FFT as a subroutine.
-
-
-'### Helper functions
-
-def odd_sized_palindrome (mid:a) (seq:n=>a) :
-  ({backward:n | mid:Unit | zforward:n}=>a) =  -- Alphabetical order matters here.
-  -- Turns sequence 12345 into 543212345.
-  for i.
-    case i of
-      {|backward=i|} -> seq.(reflect i)
-      {|mid=i|} -> mid
-      {|zforward=i|} -> seq.i
-
-
-'## Inner FFT functions
-
-data FTDirection =
-  ForwardFT
-  InverseFT
-
-def butterfly_ixs (j':halfn) (pow2:Int) : (n & n & n & n) =
-  -- Re-index at a finer frequency.
-  -- halfn must have half the size of n.
-  -- For explanation, see https://en.wikipedia.org/wiki/Butterfly_diagram
-  -- Note: with fancier index sets, this might be replacable by reshapes.
-  j = ordinal j'
-  k = ((idiv j pow2) * pow2 * 2) + mod j pow2
-  left_write_ix  = unsafeFromOrdinal n k
-  right_write_ix = unsafeFromOrdinal n (k + pow2)
-
-  left_read_ix  = unsafeFromOrdinal n j
-  right_read_ix = unsafeFromOrdinal n (j + size halfn)
-  (left_read_ix, right_read_ix, left_write_ix, right_write_ix)
-
-def power_of_2_fft (direction: FTDirection) (x: n=>Complex) : n=>Complex =
-  -- Input size must be a power of 2.
-  -- Can enforce this with tables-as-index-sets like:
-  -- (x: (log2n=>(Fin 2))=>Complex)) once that's supported.
-  dir_const = case direction of
-    ForwardFT -> -pi
-    InverseFT -> pi
-
-  log2n = intlog2 (size n)
-  halfn = idiv (size n) 2
-  
-  ans = yieldState x \xRef.
-    for i:(Fin log2n).
-      ipow2 = intpow 2 (ordinal i)
-      xRef := yieldAccum (AddMonoid Complex) \bufRef.
-        for j:(Fin halfn).  -- Executes in parallel.
-          (left_read_ix, right_read_ix,
-           left_write_ix, right_write_ix) = butterfly_ixs j ipow2
-
-          -- Read one element from the last buffer, scaled.
-          angle = dir_const * (IToF $ mod (ordinal j) ipow2) / IToF ipow2
-          v = (get xRef!right_read_ix) * (MkComplex (cos angle) (sin angle))
-
-          -- Add and subtract it to the relevant places in the new buffer.
-          bufRef!left_write_ix  += (get (xRef!left_read_ix)) + v
-          bufRef!right_write_ix += (get (xRef!left_read_ix)) - v
-
-  case direction of
-    ForwardFT -> ans
-    InverseFT -> ans / (IToF (size n))
-
-def convolve_complex (u:n=>Complex) (v:m=>Complex) :
-  ({orig_vals:n | padding:m }=>Complex) =  -- Alphabetical order matters here.
-  -- Convolve by pointwise multiplication in the Fourier domain.
-  convolved_size = (size n) + (size m) - 1
-  working_size = nextpow2 convolved_size
-  u_padded = padTo (Fin working_size) zero u
-  v_padded = padTo (Fin working_size) zero v
-  spectral_u = power_of_2_fft ForwardFT u_padded
-  spectral_v = power_of_2_fft ForwardFT v_padded
-  spectral_conv = for i. spectral_u.i * spectral_v.i
-  padded_conv = power_of_2_fft InverseFT spectral_conv
-  slice padded_conv 0 {orig_vals:n | padding:m }
-
-def convolve (u:n=>Float) (v:m=>Float) : ({orig_vals:n | padding:m }=>Float) =
-  u' = for i. MkComplex u.i 0.0
-  v' = for i. MkComplex v.i 0.0
-  ans = convolve_complex u' v'
-  for i.
-    (MkComplex real imag) = ans.i
-    real
-
-
-'## FFT Interface
-
-def fft (x: n=>Complex): n=>Complex =
-  if isPowerOf2 (size n)
-    then power_of_2_fft ForwardFT x
-    else
-      -- Bluestein's algorithm.
-      -- Converts the general FFT into a convolution,
-      -- which is then solved with calls to a power-of-2 FFT.
-      im = MkComplex 0.0 1.0
-      wks = for i.
-        i_squared = IToF $ sq $ ordinal i
-        exp $ (-im) * (MkComplex (pi * i_squared / (IToF (size n))) 0.0)
-
-      (AsList _ tailTable) = tail wks 1
-      back_and_forth = odd_sized_palindrome (head wks) tailTable
-      xq = for i. x.i * wks.i
-      back_and_forth_conj = for i. complex_conj back_and_forth.i
-      convolution = convolve_complex xq back_and_forth_conj
-      convslice = slice convolution (size n - 1) n
-      for i. wks.i * convslice.i
-
-def ifft (xs: n=>Complex): n=>Complex =
-  if isPowerOf2 (size n)
-    then power_of_2_fft InverseFT xs
-    else
-      unscaled_fft = fft (for i. complex_conj xs.i)
-      for i. (complex_conj unscaled_fft.i) / (IToF (size n))
-
-def  fft_real (x: n=>Float): n=>Complex =  fft for i. MkComplex x.i 0.0
-def ifft_real (x: n=>Float): n=>Complex = ifft for i. MkComplex x.i 0.0
-
-def fft2 (x: n=>m=>Complex): n=>m=>Complex =
-  x'      = for i. fft x.i
-  transpose for i. fft (transpose x').i
-
-def ifft2 (x: n=>m=>Complex): n=>m=>Complex =
-  x'      = for i. ifft x.i
-  transpose for i. ifft (transpose x').i
-
-def  fft2_real (x: n=>m=>Float): n=>m=>Complex =  fft2 for i j. MkComplex x.i.j 0.0
-def ifft2_real (x: n=>m=>Float): n=>m=>Complex = ifft2 for i j. MkComplex x.i.j 0.0
-
--------- Tests --------
-
-a = for i. MkComplex [10.1, -2.2, 8.3, 4.5, 9.3].i 0.0
-b = for i:(Fin 20) j:(Fin 70).
-  MkComplex (randn $ ixkey2 (newKey 0) i j) 0.0
-
-:p a ~~ (ifft $ fft a)
-> True
-:p a ~~ (fft $ ifft a)
-> True
-:p b ~~ (ifft2 $ fft2 b)
-> True
-:p b ~~ (fft2 $ ifft2 b)
-> True
diff --git a/lib/parser.dx b/lib/parser.dx
index 1190fe84e..e3d8b18d0 100644
--- a/lib/parser.dx
+++ b/lib/parser.dx
@@ -42,6 +42,12 @@ def pChar (c:Char) : Parser Unit = MkParser \(s, posRef).
   assert (c == c')
   posRef := i + 1
 
+def pString (x:String) : Parser Unit = MkParser \h.
+  (AsList n ls) = x
+  for i : (Fin n).
+      parse h (pChar ls.i)
+  ()
+
 def pEOF : Parser Unit = MkParser \(s, posRef).
   assert $ get posRef >= listLength s
 
diff --git a/lib/prelude.dx b/lib/prelude.dx
index 2879e9131..f22407ec6 100644
--- a/lib/prelude.dx
+++ b/lib/prelude.dx
@@ -1462,8 +1462,19 @@ def withFile (f:FilePath) (mode:StreamMode)
 def writeFile (f:FilePath) (s:String) : {IO} Unit =
   withFile f WriteMode \stream. fwrite stream s
 
+
+
 def readFile (f:FilePath) : {IO} String =
-  withFile f ReadMode \stream. fread stream
+  withFile f ReadMode (\stream.
+     yieldState (AsList _ []) \results.
+       iter \_.
+           str = fread stream
+           (AsList n _) = str
+           results := (get results) <> str
+           case n == 0 of
+                True -> Done ()
+                False -> Continue)
+
 
 def hasFile (f:FilePath) : {IO} Bool =
   stream = fopen f ReadMode
@@ -1997,6 +2008,55 @@ def softmax (x: n=>Float) : n=>Float =
   s = sum e
   for i. e.i / s
 
+'## Windowing and Generic Convolutions
+
+' These functions support typed windowing (possibly non-linear) and
+  linear convolution style operations. Not all windowing operations
+  are elegant under the dex type system, as the standard "valid"
+  convolution will change the size of the output in non-trivial ways.
+  We therefore focus on the "same" convolution that assumes the center
+  of the window is always within the original table, thus keeping the
+  same size. This requires the user to provide enough padding for the
+  window to work.
+
+def Pad (left:Type) (mid:Type) (right:Type) : Type =
+        {left:left | mid:mid | right:right}
+def Window (left:Type) (right:Type) : Type =
+        {left:left | mid:Unit | right:right}
+
+def window (tab: Pad left m right => n) :
+           m => (Window left right) => n =
+    for i : m.
+        for j : (Window left right).
+            k = fromOrdinal _ $ (ordinal i) + (ordinal j)
+            tab.k
+
+def pad (init:n) (tab: m=>n) : (Pad left m right) => n =
+    for i. case i of
+             {|mid=j|} -> tab.j
+             {|left=j|} -> init
+             {|right=j|} -> init
+
+' Needs to  be called with `castTable` to do the striding.
+
+def stride (tab : (m & Fin len)=>n) : m => n =
+    for i. tab.(i, 0@_)
+
+' Two dimensional versions of these functions for convenience.
+
+def window2 (tab: Pad up m down => Pad left n right => o) :
+             m => n => (Window up down) => (Window left right) => o =
+    x  = window tab
+    x' = for i j. window x.i.j
+    for i j k l. x'.i.k.j.l
+
+def pad2 (init:o) (tab: m=>n=>o) : (Pad up m down) => (Pad left n right) => o =
+         pad (for i. init) (for j. pad init tab.j)
+
+def stride2 (tab : (m & Fin vlen)=>(n & Fin hlen) => o) : m => n => o=
+    for i j. tab.(i, 0@_).(j, 0@_)
+
+
 '## Polynomials
 TODO: Move this somewhere else
 

From 58dcfb3c3fdc8a39c2a472e8debffe9c8fd22ebe Mon Sep 17 00:00:00 2001
From: Sasha Rush <srush.research@gmail.com>
Date: Wed, 21 Jul 2021 19:06:43 -0400
Subject: [PATCH 5/5] update

---
 makefile | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/makefile b/makefile
index a48a47f53..d126ef5b9 100644
--- a/makefile
+++ b/makefile
@@ -100,7 +100,7 @@ example-names = mandelbrot pi sierpinski rejection-sampler \
                 regression brownian_motion particle-swarm-optimizer \
                 ode-integrator mcmc ctc raytrace particle-filter \
                 isomorphisms ode-integrator fluidsim \
-                sgd chol fft tutorial vega-plotting kernelregression \
+                sgd chol tutorial vega-plotting kernelregression \
                 quaternions manifold-gradients schrodinger
 
 test-names = uexpr-tests adt-tests type-tests eval-tests show-tests \