import Complex
-- Output sequences of Sigma/Delta modulator (1st order non-leaky
-- accumulator) approximating constant rational.

bit False = 0    
bit True  = 1

-- Bit sequence
sd q p = f 0 where
  f s = b : f s'' where
    s' = s + q
    b  = bit $ s' >= p
    s'' = s' - b * p
    

-- N periods
sdn n q p = take (fromInteger $ n * p) $ sd q p    

-- Print
ps = concat . (map f) where
  f 0 = "."
  f _ = "X"
  
-- Build a spectrum: wave forms for values from 0/p to p/p  
spec n p = concat $ map (++ "\n") lines where
  lines = map (\q -> ps $ sdn n q p) [0..p]
s n p = putStr $ spec n p


-- Low frequency noise.  The idea is to see if the spectrum of the
-- periodic approximation signal follows the approximate linear
-- analysis, which is the inverse of the feedback filter(=integrator).

i = 0 :+ 1

-- Compute magnitude of lowest harmonic in DFT.
lfn p bs = magnitude $ (1/p') * (sum $ zipWith f bs [0..]) where
  p' = fromInteger p
  f b n = b' * exp (w * n') where
    w = i * 2 * pi / p'
    n' = fromInteger n
    b' = fromInteger b

-- Test function.  Gives the amplitude of DFT bin 1 for period = p
-- corresponding to wave form q.  This should be rougly proportional
-- to 1/p.
    
-- To get a good visual image of this, it should be plotted on a
-- log/log plot, with all p-1 representatives for each p.  Apparently
-- they do not all have the same noise amplitude.
    
tf q p = (lfn p $ sdn 1 q p)

-- Closed form expression for indices of non-zero sequence elements.
nz' q p = map n [1..q] where
  n k = ceiling $ k * p / q

-- Type conversion
nz q p = nz' (i q) (i p) where i = fromIntegral

-- Closed form expression for harmonics of S/D approximation error.
harm q p h = (sum $ map w $ nz q p) / p'' where
  p'' = fromIntegral p
  w n = exp $ (0 :+ (2 * pi * n' * h' / p')) where
        [n',h',p'] = map fromIntegral [n,h,p]