## Archive for the ‘Signal Processing’ Category

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.

### 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.

### 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.