Initial untested implementation of spectral reasoning.

diff --git a/OFC/Spectral.hs b/OFC/Spectral.hs
index 9c07fd49723faa925f3e80905b036729ed60e841..972d43ac84403d9633c4a8828b96029d4e833007 100644 (file)
--- a/OFC/Spectral.hs
+++ b/OFC/Spectral.hs
@@ -5,7 +5,8 @@ module OFC.Spectral
   , frequencyMultiple
   , Band(..)
   , negateBand
-  , isEmptyBand
+  , nullBand
+  , isNullBand
   , containsFrequency
   , containsBand
   , Spectrum
@@ -68,6 +69,22 @@ constantFrequency = Frequency 0
 frequencyMultiple :: Integer -> Frequency
 frequencyMultiple = flip Frequency $ 0
 
+multiplyFrequency :: Frequency -> Integer -> Frequency
+multiplyFrequency (Frequency m o) i = Frequency (m*i) (o*i)
+
+divideFrequency :: Frequency -> Frequency -> Integer
+divideFrequency f1 f2 = fromIntegral $ sign * absDiv f1 f2
+  where
+  absDiv a b = length . takeWhile (>= 0) . tail . iterate (subtract $ abs b) $ abs a
+  sign = if (f1 >= 0 && f2 >= 0) || (f1 < 0 && f2 < 0)
+         then 1 else (-1)
+
+halveFrequency :: Frequency -> Frequency
+halveFrequency freq@(Frequency m o) =
+  if m `mod` 2 == 0
+  then Frequency (m `div` 2) (o `div` 2)
+  else error $ "Cannot halve frequency " ++ show freq
+
 data Band =
   Band Frequency Frequency 
   deriving (Eq, Show)
@@ -78,8 +95,11 @@ instance Ord Band where
 negateBand :: Band -> Band
 negateBand (Band low high) = Band (-high + 1) (-low + 1)
 
-isEmptyBand :: Band -> Bool
-isEmptyBand (Band low high) = low == high
+nullBand :: Band
+nullBand = Band 0 0
+
+isNullBand :: Band -> Bool
+isNullBand (Band low high) = low == high
 
 restrictBandToPositive :: Band -> Band
 restrictBandToPositive (Band low high) = Band (max low 0) (max high 0)
@@ -95,14 +115,19 @@ overlaps :: Band -> Band -> Bool
 overlaps b1@(Band low1 _) b2@(Band low2 _) = 
   containsFrequency b1 low2 || containsFrequency b2 low1
 
+bandIntersection :: Band -> Band -> Band
+bandIntersection b1@(Band l1 h1) b2@(Band l2 h2) = if overlaps b1 b2 
+  then Band (max l1 l2) (min h1 h2)
+  else nullBand
+
 data Spectrum = Spectrum {
-  numSamples :: Frequency,
+  samplingRate :: Frequency,
   content :: [(Band, ContentType)]
 }
 
 constructEmptySpectrum :: Integer -> Bool -> Spectrum
 constructEmptySpectrum multiple oddSamples = Spectrum {
-  numSamples = (frequencyMultiple $ 2*multiple) + 
+  samplingRate = (frequencyMultiple $ 2*multiple) + 
                (constantFrequency $ if oddSamples then 1 else 0),
   content = [(Band (constantFrequency 0) (frequencyMultiple multiple), Zero)]
 }
@@ -110,7 +135,7 @@ constructEmptySpectrum multiple oddSamples = Spectrum {
 mergeContent :: [(Band, ContentType)] -> [(Band, ContentType)]
 mergeContent content = zip bands $ findContent <$> bands
   where
-  bands = removeOverlaps . map fst $ content
+  bands = splitOverlaps . map fst $ content
   findContent band = mconcat $ snd <$> filter (overlaps band . fst) content
 
 addContent :: [(Band, ContentType)] -> [(Band, ContentType)] -> [(Band, ContentType)]
@@ -118,13 +143,48 @@ addContent content1 content2 = mergeContent $ content1 ++ content2
 
 multiplyBands :: Band -> Band -> [Band]
 multiplyBands (Band low1 high1) (Band low2 high2) = 
-  sort . nub $ filter (not . isEmptyBand) (restrictBandToPositive <$> allResultBands)
+  sort . nub $ filter (not . isNullBand) (restrictBandToPositive <$> allResultBands)
   where
   addedAndSubtracted = [ Band (low1 + low2) (high1 + high2 - 1)
                        , Band (low1 - high2 + 1) (high1 - low2)
                        ]
   allResultBands = addedAndSubtracted ++ (negateBand <$> addedAndSubtracted)
 
+wrapContent :: Frequency -> (Band, ContentType) -> [(Band, ContentType)]
+wrapContent samplingRate (Band low high, content) = 
+  mergeContent $ undamaged:(damagedPositive ++ damagedNegative)
+  where
+  nyquist = halveFrequency samplingRate
+  numFreqs = nyquist + 1
+  negNyquist = - nyquist
+  undamaged = (bandIntersection (Band low high) (Band negNyquist numFreqs), content)
+  damagedPositive = damagedBand <$> wrapBand samplingRate (Band (max numFreqs low) (max numFreqs high))
+  damagedNegative = damagedBand <$> wrapBand samplingRate (Band (min negNyquist low) (min negNyquist high))
+  damagedBand band = (band, Damaged)
+
+wrapBand :: Frequency -> Band -> [Band]
+wrapBand samplingRate (Band low high) =
+  if samplingRate <= 0
+    then error "wrapBand takes a sampling rate greater than 0"
+    else mergeOverlaps wrappedBands
+  where
+  wrappedBands = (wrapBand' positiveBand) ++ (negateBand <$> (wrapBand' (negateBand negativeBand)))
+  positiveBand = Band (max 0 low) (max 0 high)
+  negativeBand = Band (min 0 low) (min 0 high)
+  wrapBand' (Band low high) = unwrapped : wrapped : laterBands 
+    where
+    laterBands = if localHigh == newHigh
+      then []
+      else wrapBand' (Band 0 (newHigh - samplingRate))
+    nyquist = halveFrequency samplingRate
+    numFreqs = nyquist + 1
+    offsetFrequency = samplingRate `multiplyFrequency` (low `divideFrequency` samplingRate)
+    newLow = low - offsetFrequency
+    newHigh = high - offsetFrequency
+    localHigh = min newHigh samplingRate
+    unwrapped = Band (min newLow numFreqs) (min newHigh numFreqs)
+    wrapped = Band (min (samplingRate - localHigh + 1) numFreqs) (min (samplingRate - newLow + 1) numFreqs)
+
 multiplyContent :: [(Band, ContentType)] -> [(Band, ContentType)] -> [(Band, ContentType)]
 multiplyContent bandContent1 bandContent2 = mergeContent $ do
   (band1, content1) <- bandContent1
@@ -132,8 +192,16 @@ multiplyContent bandContent1 bandContent2 = mergeContent $ do
   resultBand <- multiplyBands band1 band2
   return (resultBand, mappend content1 content2)
 
-removeOverlaps :: [Band] -> [Band]
-removeOverlaps = sort . foldl' addBand []
+mergeOverlaps :: [Band] -> [Band]
+mergeOverlaps = mergeContiguous . sort . splitOverlaps
+  where
+  mergeContiguous ((Band l1 h1):(Band l2 h2):bs) = if h1==l2
+    then (Band l1 h2):(mergeContiguous bs)
+    else (Band l1 h2):(Band l2 h2):(mergeContiguous bs)
+  mergeContiguous bs = bs
+
+splitOverlaps :: [Band] -> [Band]
+splitOverlaps = sort . foldl' addBand []
   where
   addBand bands band = nonoverlapping ++ (createBands edges) 
     where