module DiscreteWavelet.Lifting where

{- ghci -lgslcblas -lgsl -i:../math:../gslmatrix/ ../gslmatrix/gslaux.o ShiftedSignal.hs DiscreteWavelet.Lifting -}

import LinearAlgebra.Matrix(Matrix,matrix,matrix_x_matrix)

import Numerics.SpecialFunctions (factorial, recipBinomial, binomial)
import qualified ShiftedSignal as SS

{- * Lifting steps for weighting of channels -}

skewAdditionMatrix :: Double -> Matrix
skewAdditionMatrix x = matrix [[x,1],[1,0]]

additionWeightingGeneric :: Double -> Double -> [Double]
additionWeightingGeneric a x =
   let b = (x-1)/a
   in  reverse [a, b, -a/x, -b*x]

additionWeightingRational :: Double -> [Double]
additionWeightingRational x =
   reverse [1, x-1, -1/x, x*(1-x)]

-- does not work for x>1
additionWeightingOptimal :: Double -> [Double]
additionWeightingOptimal x =
   let steps y =
           let q = sqrt((1-y)/(2+y))
               p = sqrt((1-y)*(2+y))
               s = sqrt y
           in  [-p*s, q/s, p/s, -q*s]
   in if x<=1
        then reverse (steps x)
        else map negate (steps (recip x))

testAdditions :: [Double] -> Matrix
testAdditions xs =
   foldl1 matrix_x_matrix (map skewAdditionMatrix xs)



{- * Lifting steps for CDF lowpasses, allowing also fractional orders -}

{- |
  If the lifting sequence is reversed it can be used to construct
  the lowpass or the highpass independent from the other filter.
  This function generates a list of modified liftings steps
  where the target filter of convolve-accumulate operation is also weighted.
  This sequence was designed for the forward application order
  but when applied in reversed order it does not work. -}
cdfHalfLiftingSeqEven :: Fractional a => a -> [(a, SS.T a)]
cdfHalfLiftingSeqEven n =
   tail $
   zip (zipWith (*)
           (iterate (subtract 2) (2*n))
           (iterate (2+) 0))
       (zipWith SS.translate (cycle [-1,0])
           (map (SS.fromList . replicate 2)
                (iterate (subtract 2) (n+1))))

{- |
  Do some modified lifting steps where the target filter is weighted. -}
noUpsampleWeightLifting2 :: Num a => [(a, SS.T a)] -> (SS.T a, SS.T a) -> [(SS.T a, SS.T a)]
noUpsampleWeightLifting2 lseq start =
   scanl (\(x, y) (k, s) -> (y, SS.superpose (SS.amplify k y) (SS.convolve s x)))
         start lseq

{- |
  Generate lifting filters for CDF low-pass of even order.
  They are designed in a way that also allows fractional orders
  but the sequence of lifted filters does not converge. -}
cdfLiftingSeqEvenFrac :: Double -> [SS.T Double]
cdfLiftingSeqEvenFrac n =
   zipWith SS.translate
           (cycle [0,-1])
           (map (\c -> SS.fromList [c,c])
                (cdfLiftingSeqEvenFracCoeffs n))

cdfLiftingSeqEvenFracCoeffs :: Double -> [Double]
cdfLiftingSeqEvenFracCoeffs n =
   let am  = recipBinomial (n/2-0.5) 0.5 / 2
       modds  = iterate (subtract 2) (n-1)
       mevens = iterate (subtract 1) (n/2-1)
       as = scanl (\ai m -> recip ((n^(2::Int) - 4 * m^(2::Int)) * ai))
                  am mevens
   in  zipWith (*) modds as

{- factorial (n/2-1) / (4 * factorial (n/2-0.5) / factorial (-0.5))
      == factorial (n/2-1) * factorial (-0.5) / (4 * factorial (n/2-0.5))
      == factorial (n/2-1) * factorial   0.5  / (2 * factorial (n/2-0.5)) -}

{- The lowpass as it evolves from lifting steps -}
cdfLiftingEvenFracLowpass :: Double -> [SS.T Double]
cdfLiftingEvenFracLowpass n =
   map (SS.interleave . (\(x,y) -> [y,x]))
       (scanl (flip SS.noUpsampleLiftingStep2)
              (SS.fromScalar 1, SS.fromScalar 0)
              (cdfLiftingSeqEvenFrac n))

{-
  mapM (\x -> plotList [] (ShiftedSignal.signal (ShiftedSignal.multiRefine 6 2 x (ShiftedSignal.fromScalar 1)))) (take 10 $ ShiftedSignal.cdfLiftingEvenFracLowpass 4.3)
-}