-----------------------------------------------------------------------------
-- |
-- Module      :  DSP.Multirate.Polyphase
-- Copyright   :  (c) Matthew Donadio 2003
-- License     :  GPL
--
-- Maintainer  :  m.p.donadio@ieee.org
-- Stability   :  experimental
-- Portability :  portable
--
-- Polyphase interpolators and decimators
--
-- Reference: C&R
--
-----------------------------------------------------------------------------

module DSP.Multirate.Polyphase (poly_interp) where

import Data.List (transpose)
import Data.Array

import DSP.Filter.FIR.FIR

-- mkpoly turns a single filter into a list of l subfilters

mkpoly :: Num a => Array Int a -> Int -> Int -> Array Int a
mkpoly :: forall a. Num a => Array Int a -> Int -> Int -> Array Int a
mkpoly Array Int a
h Int
l Int
k = forall i e. Ix i => (i, i) -> [e] -> Array i e
listArray (Int
0,Int
m) [ Array Int a
hforall i e. Ix i => Array i e -> i -> e
!(Int
kforall a. Num a => a -> a -> a
+Int
nforall a. Num a => a -> a -> a
*Int
l) | Int
n <- [Int
0..Int
m] ]
    where m :: Int
m = ((forall a b. (a, b) -> b
snd forall a b. (a -> b) -> a -> b
$ forall i e. Array i e -> (i, i)
bounds Array Int a
h) forall a. Num a => a -> a -> a
+ Int
1) forall a. Integral a => a -> a -> a
`div` Int
l forall a. Num a => a -> a -> a
- Int
1

-- | Polyphase interpolator

poly_interp :: (Num a, Eq a) => Int -- ^ L
	    -> Array Int a -- ^ h[n]
	    -> [a] -- ^ x[n]
	    -> [a] -- ^ y[n]

poly_interp :: forall a. (Num a, Eq a) => Int -> Array Int a -> [a] -> [a]
poly_interp Int
l Array Int a
h [a]
x = forall (t :: * -> *) a. Foldable t => t [a] -> [a]
concat forall a b. (a -> b) -> a -> b
$ forall a. [[a]] -> [[a]]
transpose [[a]]
y
    where g :: [[a] -> [a]]
g = forall a b. (a -> b) -> [a] -> [b]
map (forall a. (Num a, Eq a) => Array Int a -> [a] -> [a]
fir forall b c a. (b -> c) -> (a -> b) -> a -> c
. forall a. Num a => Array Int a -> Int -> Int -> Array Int a
mkpoly Array Int a
h Int
l) [Int
0..(Int
lforall a. Num a => a -> a -> a
-Int
1)]
	  y :: [[a]]
y = forall a b. (a -> b) -> [a] -> [b]
map (forall a b. (a -> b) -> a -> b
$ [a]
x) [[a] -> [a]]
g

{-

gZipWith :: Eq a => (a -> a -> a) -> [[a]] -> [a]
gZipWith f xs | any (== []) xs = []
	      | otherwise = foldl1 f (map head xs) : gZipWith f (map tail xs)

poly_decim :: Num a => Int -> Array Int a -> [a] -> [a]

poly_decim l h x = gZipWith (+) g
    where g = map (fir . mkpoly h l) [0..(l-1)]

Test

> h :: Array Int Double
> h = listArray (0,15) [1..16]

> x :: [Double]
> x =  [ 1, 0, 0, 0 ]

-}