--- /dev/null
+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 _ = []