From: Francis Russell <francis@unchartedbackwaters.co.uk>
Date: Fri, 21 Dec 2012 14:48:04 +0000 (+0000)
Subject: Add initial work on spectrum analysis.
X-Git-Url: https://git.unchartedbackwaters.co.uk/w/?a=commitdiff_plain;h=29609ec4a2c9997e23cc2009b2be627083e1d7de;p=francis%2Fofc.git

Add initial work on spectrum analysis.
---

diff --git a/OFC/Spectral.hs b/OFC/Spectral.hs
new file mode 100644
index 0000000..9c07fd4
--- /dev/null
+++ b/OFC/Spectral.hs
@@ -0,0 +1,145 @@
+module OFC.Spectral
+  ( ContentType(..)
+  , Frequency
+  , constantFrequency
+  , frequencyMultiple
+  , Band(..)
+  , negateBand
+  , isEmptyBand
+  , containsFrequency
+  , containsBand
+  , Spectrum
+  , constructEmptySpectrum
+  ) where
+
+import OFC.Common(PrettyPrintable(..))
+
+import Control.Applicative ((<$>))
+import Data.List (foldl', nub, sort, partition)
+import Data.Monoid (Monoid(..))
+import Text.PrettyPrint (text)
+
+data ContentType 
+  = Zero
+  | Present
+  | Damaged
+  deriving (Eq, Show)
+
+instance Monoid ContentType where
+  mempty = Zero
+  mappend Zero s = s
+  mappend s Zero = s
+  mappend Damaged _ = Damaged
+  mappend _ Damaged = Damaged
+  mappend Present Present = Present
+
+type Multiple = Integer
+type Offset = Integer
+
+data Frequency 
+  = Frequency Multiple Offset
+  deriving (Eq, Show)
+
+instance Num Frequency where
+  (+) (Frequency m1 o1) (Frequency m2 o2) = Frequency (m1+m2) (o1+o2)
+  (*) _ _ = error "Multiplication undefined for Frequency"
+  negate (Frequency m o) = Frequency (-m) (-o)
+  abs freq = if freq < fromInteger 0
+    then negate freq 
+    else freq
+  signum _ =  error "signum undefined for Frequency"
+  fromInteger = constantFrequency
+
+instance Ord Frequency where
+  compare (Frequency m1 o1) (Frequency m2 o2) = compare (m1, o1) (m2, o2)
+
+instance PrettyPrintable Frequency where
+  toDoc (Frequency m o) = if m == 0 
+    then text $ show o 
+    else text $ show m ++ "N" ++ constantComponent o where
+      constantComponent c
+        | c < 0  = " - " ++ (show $ abs c)
+        | c > 0  = " + " ++ show c
+        | otherwise = ""
+
+constantFrequency :: Integer -> Frequency
+constantFrequency = Frequency 0
+
+frequencyMultiple :: Integer -> Frequency
+frequencyMultiple = flip Frequency $ 0
+
+data Band =
+  Band Frequency Frequency 
+  deriving (Eq, Show)
+
+instance Ord Band where
+  compare (Band l1 h1) (Band l2 h2) = compare (l1, h1) (l2, h2)
+
+negateBand :: Band -> Band
+negateBand (Band low high) = Band (-high + 1) (-low + 1)
+
+isEmptyBand :: Band -> Bool
+isEmptyBand (Band low high) = low == high
+
+restrictBandToPositive :: Band -> Band
+restrictBandToPositive (Band low high) = Band (max low 0) (max high 0)
+
+containsFrequency :: Band -> Frequency -> Bool
+containsFrequency (Band low high) freq = freq >= low && freq < high
+
+containsBand :: Band -> Band -> Bool
+containsBand (Band low high) (Band sLow sHigh) = 
+  sLow >= low && sHigh <= high
+
+overlaps :: Band -> Band -> Bool
+overlaps b1@(Band low1 _) b2@(Band low2 _) = 
+  containsFrequency b1 low2 || containsFrequency b2 low1
+
+data Spectrum = Spectrum {
+  numSamples :: Frequency,
+  content :: [(Band, ContentType)]
+}
+
+constructEmptySpectrum :: Integer -> Bool -> Spectrum
+constructEmptySpectrum multiple oddSamples = Spectrum {
+  numSamples = (frequencyMultiple $ 2*multiple) + 
+               (constantFrequency $ if oddSamples then 1 else 0),
+  content = [(Band (constantFrequency 0) (frequencyMultiple multiple), Zero)]
+}
+
+mergeContent :: [(Band, ContentType)] -> [(Band, ContentType)]
+mergeContent content = zip bands $ findContent <$> bands
+  where
+  bands = removeOverlaps . map fst $ content
+  findContent band = mconcat $ snd <$> filter (overlaps band . fst) content
+
+addContent :: [(Band, ContentType)] -> [(Band, ContentType)] -> [(Band, ContentType)]
+addContent content1 content2 = mergeContent $ content1 ++ content2
+
+multiplyBands :: Band -> Band -> [Band]
+multiplyBands (Band low1 high1) (Band low2 high2) = 
+  sort . nub $ filter (not . isEmptyBand) (restrictBandToPositive <$> allResultBands)
+  where
+  addedAndSubtracted = [ Band (low1 + low2) (high1 + high2 - 1)
+                       , Band (low1 - high2 + 1) (high1 - low2)
+                       ]
+  allResultBands = addedAndSubtracted ++ (negateBand <$> addedAndSubtracted)
+
+multiplyContent :: [(Band, ContentType)] -> [(Band, ContentType)] -> [(Band, ContentType)]
+multiplyContent bandContent1 bandContent2 = mergeContent $ do
+  (band1, content1) <- bandContent1
+  (band2, content2) <- bandContent2
+  resultBand <- multiplyBands band1 band2
+  return (resultBand, mappend content1 content2)
+
+removeOverlaps :: [Band] -> [Band]
+removeOverlaps = sort . foldl' addBand []
+  where
+  addBand bands band = nonoverlapping ++ (createBands edges) 
+    where
+    (overlapping, nonoverlapping) = partition (overlaps band) bands
+    edges = getEdges $ band : overlapping
+    getEdges = sort . nub . concatMap getFrequencyEdges 
+    getFrequencyEdges (Band low high) = [low, high]
+    createBands (f1:f2:fs) = (Band f1 f2) : createBands (f2:fs)
+    createBands _ = []
diff --git a/ofc.cabal b/ofc.cabal
index 574813f..9e91d80 100644
--- a/ofc.cabal
+++ b/ofc.cabal
@@ -15,6 +15,7 @@ Library
                      OFC.TopLevel
                      OFC.SecondLevel
                      OFC.TargetMapping
+                     OFC.Spectral
 
 Executable ofc
   Main-is:           Main.hs