-----------------------------------------------------------------------------
-- |
-- Module      :  DSP.Convolution
-- Copyright   :  (c) Matthew Donadio 2003
-- License     :  GPL
--
-- Maintainer  :  m.p.donadio@ieee.org
-- Stability   :  experimental
-- Portability :  portable
--
-- Module to perform the linear convolution of two sequences
--
-----------------------------------------------------------------------------

module DSP.Convolution (conv, test) where

import Data.Array

-- * Functions

-- | @conv@ convolves two finite sequences

conv :: (Ix a, Integral a, Num b) => Array a b -> Array a b -> Array a b
conv :: forall a b.
(Ix a, Integral a, Num b) =>
Array a b -> Array a b -> Array a b
conv Array a b
x1 Array a b
x2 = Array a b
x3
    where m1 :: a
m1 = forall a b. (a, b) -> b
snd forall a b. (a -> b) -> a -> b
$ forall i e. Array i e -> (i, i)
bounds Array a b
x1
          m2 :: a
m2 = forall a b. (a, b) -> b
snd forall a b. (a -> b) -> a -> b
$ forall i e. Array i e -> (i, i)
bounds Array a b
x2
	  m3 :: a
m3 = a
m1 forall a. Num a => a -> a -> a
+ a
m2
	  x3 :: Array a b
x3 = forall i e. Ix i => (i, i) -> [e] -> Array i e
listArray (a
0,a
m3) [
                    forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
sum [ Array a b
x1forall i e. Ix i => Array i e -> i -> e
!a
k forall a. Num a => a -> a -> a
* Array a b
x2forall i e. Ix i => Array i e -> i -> e
!(a
nforall a. Num a => a -> a -> a
-a
k) | a
k <- [forall a. Ord a => a -> a -> a
max a
0 (a
nforall a. Num a => a -> a -> a
-a
m2)..forall a. Ord a => a -> a -> a
min a
n a
m1] ]
                       | a
n <- [a
0..a
m3] ]

-- Test vectors.  Linear convolution is also equivalent to polynomial
-- multiplication.

h1, h2, h3 :: Array Int Integer
h1 :: Array Int Integer
h1 = forall i e. Ix i => (i, i) -> [e] -> Array i e
listArray (Int
0,Int
3) [ Integer
1, Integer
2, Integer
3, Integer
4 ]
h2 :: Array Int Integer
h2 = forall i e. Ix i => (i, i) -> [e] -> Array i e
listArray (Int
0,Int
4) [ Integer
1, Integer
2, Integer
3, Integer
4, Integer
5 ]
h3 :: Array Int Integer
h3 = forall i e. Ix i => (i, i) -> [e] -> Array i e
listArray (Int
0,Int
7) [ Integer
1, Integer
4, Integer
10, Integer
20, Integer
30, Integer
34, Integer
31, Integer
20 ]

test :: Bool
test :: Bool
test  =  forall a b.
(Ix a, Integral a, Num b) =>
Array a b -> Array a b -> Array a b
conv Array Int Integer
h1 Array Int Integer
h2 forall a. Eq a => a -> a -> Bool
== Array Int Integer
h3