Signal Processing | Rod Carvalho

Archive for the ‘Signal Processing’ Category

Circular convolution in Haskell

September 5, 2011

Given two finite sequences, $x$ and $y$ , of length $n$ , their $n$ -point circular convolution can be written in matrix form [1] as follows

$\left[\begin{array}{c} z_0\\ z_1\\ z_2\\ \vdots\\ z_{n-2}\\ z_{n-1}\end{array}\right] = \left[\begin{array}{cccccc} y_0 & y_{n-1} & y_{n-2} & \ldots & y_2 & y_1\\ y_1 & y_0 & y_{n-1} & \ldots & y_3 & y_2\\ y_2 & y_1 & y_0 & \ldots & y_4 & y_3\\ \vdots & \vdots & \vdots & \ddots & \vdots & \vdots\\ y_{n-2} & y_{n-3} & y_{n-4} & \ldots & y_0 & y_{n-1}\\ y_{n-1} & y_{n-2} & y_{n-3} & \ldots & y_1 & y_0\end{array}\right] \left[\begin{array}{c} x_0\\ x_1\\ x_2\\ \vdots\\ x_{n-2}\\ x_{n-1}\end{array}\right]$

where $z$ is also a sequence of length $n$ . Note that the $n \times n$ matrix above is a circulant matrix, as each row is a circularly right-shifted version of the previous one. Each $z_k$ is obtained by taking the inner product of the $k$ -th row of the matrix with the $x$ vector.

In Digital Signal Processing (DSP), we often think of finite sequences as vectors, which allows us to compute convolutions using matrix operations. Taking a more Computer Science-ish approach, one can think of finite sequences as finite lists, and compute convolutions using list operations. The matrix above can be thought of as a list of lists (where each row is a list). We know how to multiply a matrix by a vector, but how do we “multiply” a list of lists by a list?

I started teaching myself Haskell this Summer, and I am (obviously) looking for interesting problems at which to throw some functional firepower. Solutions looking for a problem. But, first, let us “begin at the beginning”…

__________

Reverse

Take another look at the matrix above and note that its last row, when viewed as a list, is the reversal of list $y$ , which means that we need a function that reverses lists. Fortunately, we do not need to implement such a function, as there is the function reverse in Prelude that uses foldl and flip:

reverse :: [a] -> [a]
reverse = foldl (flip (:)) []

__________

Circular shifts

Yet once again, let us look at the matrix at the top of this post. Note that the 2nd row is a circularly right-shifted version of the 1st one, and that the 3rd row is a circularly right-shifted version of the 2nd one, et cetera. Also, note that the $(n-1)$ -th row is a circularly left-shifted version of the $n$ -th row, and so and so on. Thus, we need to implement functions that perform circular shifts.

The following function circularly right-shifts a list:

circShiftR :: [a] -> [a]
circShiftR [] = []
circShiftR x = last x : init x

This function merely takes the last element of the list and moves it to the beginning. We can circularly left-shift a list using the function:

circShiftL :: [a] -> [a]
circShiftL [] = []
circShiftL xs = (tail xs) ++ [head xs]

which moves the first element of a list to the end. Let us now test these functions on some simple lists using the GHCi interpreter:

*Main> circShiftR [0..9]
[9,0,1,2,3,4,5,6,7,8]
*Main> circShiftR (circShiftR [0..9])
[8,9,0,1,2,3,4,5,6,7]
*Main> (circShiftR . circShiftR) [0..9]
[8,9,0,1,2,3,4,5,6,7]
*Main> circShiftL [0..9]
[1,2,3,4,5,6,7,8,9,0]
*Main> (circShiftL . circShiftL) [0..9]
[2,3,4,5,6,7,8,9,0,1]

Note that in the 3rd and 5th inputs we composed the circular shifting functions to perform double right- and left-shifts.

__________

Iterate

Since each row of the matrix is a right-shifted version of the previous one or a left-shifted version of the next one, we can generate the matrix from its first or last rows. Putting it in terms of lists, we can generate a list of lists that represents the matrix from a single list using the iterate function in Prelude:

iterate :: (a -> a) -> a -> [a]
iterate f x =  x : iterate f (f x)

and the circular shifting functions we implemented before. To make things concrete, we can generate all four circular right-shifts of a list of length 4, as follows:

*Main> take 4 (iterate circShiftR [0..3])
[[0,1,2,3],[3,0,1,2],[2,3,0,1],[1,2,3,0]]

A word of caution is in order. Note that we take the first four lists, as

iterate f x

does generate an infinite list containing $x$ , $f (x)$ , $(f \circ f ) (x)$ , $(f \circ f \circ f ) (x)$ , et cetera. We do not need an infinite list of lists.

We now know how to generate a list of lists representing the matrix.

__________

Inner product

To obtain the circular convolution of $x$ and $y$ , we multiply a matrix by a (column) vector, which consists of computing $n$ inner products. However, in this post we are thinking in terms of lists. We thus need to implement a function that computes the inner product of two (finite) lists, as follows:

innerProd :: Num a => [a] -> [a] -> a
innerProd xs ys = sum(zipWith (*) xs ys)

where we first use zipWith to obtain a list of the $x_i y_i$ products, and then use sum to compute the sum $\sum_i x_i y_i$ . Let us test the inner product function:

*Main> let xs = [1..4]
*Main> let ys = [0..3]
*Main> zipWith (*) xs ys
[0,2,6,12]
*Main> sum(zipWith (*) xs ys)
20

__________

Circular convolution

We finally have all we need to compute the circular convolution.

For example, suppose we want to compute the circular convolution of two 4-point sequences, say, $x = (1,1,1,1)$ and $y = (0,1,2,3)$ . Using matrix-vector multiplication, we obtain:

$\left[\begin{array}{cccc} 0 & 3 & 2 & 1\\ 1 & 0 & 3 & 2\\ 2 & 1 & 0 & 3\\ 3 & 2 & 1 & 0\end{array}\right] \left[\begin{array}{c} 1\\ 1\\ 1\\ 1\end{array}\right] = \left[\begin{array}{c} 6\\ 6\\ 6\\ 6\end{array}\right]$

Alternatively, we can compute the circular convolution by representing sequences by lists and using the functions we have implemented before:

*Main> -- define lists
*Main> let xs = [1,1,1,1]
*Main> let ys = [0,1,2,3]
*Main> -- reverse and right-shift ys
*Main> let ys' = (circShiftR . reverse) ys
*Main> ys'
[0,3,2,1]
*Main> -- compute list of lists representing the matrix
*Main> let yss = take 4 (iterate circShiftR ys')
*Main> yss
[[0,3,2,1],[1,0,3,2],[2,1,0,3],[3,2,1,0]]
*Main> -- compute inner product of xs with each list in yss
*Main> map (innerProd xs) yss
[6,6,6,6]

Note that I wrote comments on the GHCi command line in order to improve readability. We are now ready to implement the function that computes the circular convolution:

circConv :: Num a => [a] -> [a] -> [a]
circConv xs ys = map (innerProd xs) yss
		 where
	     	   n   = length xs
	           ys' = (circShiftR . reverse) ys
		   yss = take n (iterate circShiftR ys')

How does it work? We first take sequence $y$ , reverse and right-shift it (to obtain the 1st row of the matrix), generate the remaining rows by right-shifting the previous ones, and finally compute the inner product of each row with sequence $x$ using map. Note that

innerProd xs

is a partial application of the inner product function, as we provide only the first argument. We could define a function $g$ as follows:

*Main> let g = innerProd [1,1,1,1]
*Main> :type g
g :: [Integer] -> Integer
*Main> g [0..3]
6

Note that $g$ takes a list of integers and returns its inner product with the list of four ones. Putting it all together in a .hs file, we have:

-- circular right-shift
circShiftR :: [a] -> [a]
circShiftR [] = []
circShiftR x = last x : init x

-- inner product of two lists
innerProd :: Num a => [a] -> [a] -> a
innerProd xs ys = sum(zipWith (*) xs ys)

-- circular convolution of two lists
circConv :: Num a => [a] -> [a] -> [a]
circConv xs ys = map (innerProd xs) yss
		 where
	     	    n   = length xs
	            ys' = (circShiftR . reverse) ys
		    yss = take n (iterate circShiftR ys')

We use several higher-order functions in this implementation: map, iterate, zipWith. We also use function composition and partial function application. A lot of beautiful ideas in only a dozen lines of code!

I am a Haskell neophyte, and I make no claims that the code above cannot be improved. If you have any constructive suggestions, I would be happy to hear them. Please use the HTML tags <pre> and </pre> to post code in the comments.

__________

Testing

Last but not least, let us carry out a couple of tests using examples from Mitra’s book [1]:

*Main> -- problem 5.2 b)
*Main> let xs = [-1,5,3,0,3]
*Main> let hs = [-2,0,5,3,-2]
*Main> circConv xs hs
[1,-1,-2,16,26]
*Main> -- problem 5.45 b)
*Main> let gs = [2,-1,3,0]
*Main> let hs = [-2,4,2,-1]
*Main> circConv gs hs
[3,7,-6,8]

It appears to be working! Please let me know if you find any bugs.

__________

References

[1] Sanjit K. Mitra, Digital Signal Processing: a computer-based approach, 3rd edition, McGraw-Hill, 2005.

Tags:Circular Convolution, Digital Signal Processing, Discrete-Time Signal Processing, Functional Programming, Haskell, Higher-Order Functions, List Processing
Posted in Haskell, Signal Processing | 2 Comments »

A convoluted convolution

April 2, 2011

Consider the signal $x(t) = 1/t$ . What is the convolution of $x$ with itself?

__________

My solution:

From the definition of convolution [1], we have

$(x * x) (t) = \displaystyle \int_{-\infty}^{+\infty}\frac{d \tau}{\tau (t - \tau)} = \displaystyle \frac{1}{t} \int_{-\infty}^{+\infty} \left(\frac{1}{\tau} + \frac{1}{t - \tau}\right) d \tau$ .

Unfortunately, the latter integral is plagued by convergence problems. To make matters worse, $x(t)$ is not even defined at $t = 0$ . Fortunately, the former integral is a familiar one: it’s a scaled Hilbert transform of $x$ . It is well-known [1] that the Fourier transform of $1 / \pi t$ (the impulse response of a Hilbert transformer) is the following

$\mathcal{F} \left(\displaystyle\frac{1}{\pi t}\right) = - j \,\text{sgn} (w) =\begin{cases} -j, & \omega > 0\\ +j, & \omega < 0\end{cases}$

where $\text{sgn}(\cdot)$ is the signum function. Thus, we have that the Fourier transform of $x$ is

$\mathcal{F} (x) = -j \pi \, \text{sgn} (w) = \begin{cases} -j \pi, & \omega > 0\\ +j \pi, & \omega < 0\end{cases}$

Convolution in the time domain corresponds to multiplication in the frequency domain. Hence, the Fourier transform of the convolution is

$\mathcal{F} (x * x) = \left(\mathcal{F} (x)\right)^2 = \left(-j \pi \, \text{sgn} (w)\right)^2 = - \pi^2$

and, taking the inverse Fourier transform, we finally obtain

$(x * x) (t) = \mathcal{F}^{-1} (- \pi^2) = -\pi^2 \delta (t)$

where $\delta (t)$ is the Dirac delta. This is a prime example of a problem that Signal Processing professors in engineering departments tend to love, but that drives (some) mathematicians up the wall.

__________

References

[1] Alan V. Oppenheim, Alan S. Willsky, S. Hamid Nawab, Signals & Systems, 2nd edition, Prentice-Hall, 1996.

Tags:Hilbert Transform, Puzzles
Posted in Puzzles, Signal Processing | 2 Comments »

Interpolation and energy amplification

March 25, 2011

Consider the following interpolator [1]

where $x_u = (\uparrow L) (x)$ is the upsampling of sequence $x$ by an integer factor of $L$ , and where $H(z)$ is the transfer function of an ideal discrete-time LTI lowpass filter with gain $L$ and cutoff frequency $\pi/L$ . What is the relation between the energies of signals $x$ and $y$ ?

__________

My solution:

It can be easily shown [1] that the interpolator’s input-output relationship is $Y (z) = H(z) X (z^L)$ . Hence, evaluating the Z-transforms on the unit circle, we obtain

$Y (e^{j \omega}) = H(e^{j \omega}) X (e^{j \omega L})$ .

The energy of the interpolator’s output signal is

$\mathcal{E}_y := \displaystyle \sum_{n \in \mathbb{Z}} | y(n) |^2 = \frac{1}{2 \pi} \int_{-\pi}^{\pi} |Y (e^{j \omega})|^2 d \omega$

where we used Parseval’s theorem for discrete-time signals. Thus,

$\mathcal{E}_y = \displaystyle\frac{1}{2 \pi} \int_{-\pi}^{\pi} |H(e^{j \omega})|^2 |X (e^{j \omega L})|^2 d \omega$ .

$H(z)$ is an ideal lowpass filter with the following frequency response

$H(e^{j \omega}) = \begin{cases} L, & |\omega| \leq \pi/L\\ 0, & \pi/ L< |\omega| \leq \pi\end{cases}$

and the following sinc impulse response

$h(n) = \displaystyle\frac{1}{2 \pi} \int_{-\pi}^{\pi} H(e^{j \omega}) e^{j \omega n} d \omega = \displaystyle \frac{\sin(\frac{n \pi}{L})}{\frac{n \pi}{L}} =: \text{sinc} \left(\frac{n \pi}{L}\right)$ .

Please do note that this is a non-causal filter with an impulse response of infinite length. We can then write the interpolator’s output signal’s energy as follows

$\mathcal{E}_y = \displaystyle\frac{1}{2 \pi} \int_{-\pi}^{\pi} |H(e^{j \omega})|^2 |X (e^{j \omega L})|^2 d \omega = \displaystyle\frac{L^2}{2 \pi} \int_{-\pi/L}^{\pi/L} |X (e^{j \omega L})|^2 d \omega$ .

Performing a change of variable, $\theta = L \omega$ , we obtain

$\mathcal{E}_y = \displaystyle\frac{L^2}{2 \pi} \int_{-\pi/L}^{\pi/L} |X (e^{j \omega L})|^2 d \omega = \displaystyle\frac{L}{2 \pi} \int_{-\pi}^{\pi} |X (e^{j \theta})|^2 d \theta = L \mathcal{E}_x$

where

$\mathcal{E}_x := \displaystyle \sum_{n \in \mathbb{Z}} | x(n) |^2 = \frac{1}{2 \pi} \int_{-\pi}^{\pi} |X (e^{j \omega})|^2 d \omega$

is the energy of the input signal. Since $\mathcal{E}_y = L \mathcal{E}_x$ , we can conclude that the interpolator amplifies the input signal’s energy by a factor of $L$ , which is somewhat intuitive.

We considered the case where signals are bi-infinite sequences and $H(z)$ is a non-causal ideal lowpass filter. It would be interesting to consider the more realistic case where the signals are finite sequences and $H(z)$ is a causal FIR filter.

__________

References

[1] Alan V. Oppenheim, Ronald W. Schafer, Discrete-Time Signal Processing, 2nd edition, Prentice-Hall, 1999.

Tags:Discrete-Time Signal Processing, Interpolation, Multirate Signal Processing, SE-S2010, Signal Processing
Posted in Numerical Methods, Puzzles, Signal Processing | 2 Comments »

S	M	T	W	T	F	S
« Dec
			1	2	3	4
5	6	7	8	9	10	11
12	13	14	15	16	17	18
19	20	21	22	23	24	25
26	27	28	29	30	31

Rod Carvalho

Archive for the ‘Signal Processing’ Category

A convoluted convolution

Interpolation and energy amplification

Pages

Categories

Top Posts

Recent Posts

Recent Comments

Tag Cloud

Archives

Blog Stats

Artists

Blogs

Haskell

Math Blogs

MATLAB

Python

Robotics

Software

Symbolic Computation

Powered by FeedBurner

Meta

Follow “Rod Carvalho”

Rod Carvalho

Archive for the ‘Signal Processing’ Category

Circular convolution in Haskell

A convoluted convolution

Interpolation and energy amplification

Pages

Categories

Top Posts

Recent Posts

Recent Comments

Tag Cloud

Archives

Blog Stats

Artists

Blogs

Haskell

Math Blogs

MATLAB

Python

Robotics

Software

Symbolic Computation

Category Cloud

Powered by FeedBurner

Meta

Follow “Rod Carvalho”