Description
Basic Ideas
Consider the following figure.
0 0.2 0.4 0.6 0.8 1
time
-1.5
-1
-0.5
0
0.5
1
1.5
amplitude
corrupted
original
1
This is a simulation of a noisy signal. The idea is that a broadcasting station (e.g., a radio transmitter) has
sent out a smooth, time-varying signal s(t) (shown with the black curve). However, because of noise in the
environment (e.g., scattering from trees, interfering signals, etc.), what we have received is a corrupted version
of the signal scor(t) (the blue line). The noisy signal follows the general shape of the original signal, but also
jumps up and down rapidly in time. Our job is to reconstruct as well as possible the original signal from the
noisy signal by filtering out the noise.
Our first step is to convert the problem into the language of linear algebra. To do this, we will discretize time.
That is, rather than working with the values of s(t) and scor(t) for all real numbers t, we’ll sample s and scor
at a finite set of time values {t1, t2, . . . , tn}. That allows us to represent s and scor as vectors in R
n.
x =
s(t1)
s(t2)
.
.
.
s(tn)
and xcor =
scor(t1)
scor(t2)
.
.
.
scor(tn)
.
Our goal, then, is to use the corrupted signal xcor to produce a smoothed signal y that is hopefully close to the
original signal x (y ≈ x). One way to do this would be to use orthogonal projection onto a suitable basis. This
time, however, we’ll use a different technique, called regularization.
Let δ > 0 be a positive number (called the regularization parameter ), and, for every y in R
n, define
f(y) = ky − xcork
2 + δ
2
nX−1
i=1
(yi+1 − yi)
2
.
The first term in f(y) measures how close y is to the received signal xcor; it is small for y ≈ xcor and large for
y 6≈ xcor. The second term measures how “rough” y is. If each entry yi+1 is not too different from the previous
yi
, the sum of squares will be small. But if there are major changes from some yi to the next yi+1, these will
cause the second term of f(y) to be large.
Now, consider choosing y to make f(y) as small as possible. By our previous analysis, we should expect the best
y to strike a balance between being close to the received signal xcor and being not-too-rough. The parameter δ
quantifies the trade-off between these two goals. If δ is very close to 0, we won’t care too much about roughness
and so y will be very close to xcor. If δ is very very large, we will focus almost entirely on roughness, and make
each value of yi+1 approximately equal to the previous:
yn ≈ yn−1 ≈ yn−2 ≈ · · · ≈ y2 ≈ y1.
By tuning δ correctly, we might hope to balance between these two extremes and find a y that’s not too rough
but also not too far from xcor.
Now, before we jump into the code, we only need to recognize the problem of minimizing f(y) as a least squares
problem. To do this, we define the (n − 1) × n matrix D as
D =
−1 1 0 · · · 0
0 −1 1 · · · 0
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
0 · · · 0 −1 1
so that Dy =
y2 − y1
y3 − y2
.
.
.
yn − yn−1
.
Then
f(y) = ky − xcork
2 + δ
2
kDyk
2
.
Now we will apply one final standard trick from linear algebra. (You don’t need to understand why this works
for the lab, but thinking about it after the lab is highly recommended. In working out why the following is true,
you will be practicing a number of important skills from this class.) If we define the (2n − 1) × n matrix A and
2
the (2n − 1)-vector b as
A =
In×n
δD
=
1 0 0 · · · 0
0 1 0 · · · 0
0 0 1 · · · 0
.
.
.
.
.
.
.
.
.
.
.
.
0 0 0 · · · 1
−δ δ 0 · · · 0
0 −δ δ · · · 0
.
.
.
.
.
.
.
.
.
.
.
.
0 · · · 0 −δ δ
and b =
xcor
0n−1
=
scor(t1)
scor(t2)
scor(t3)
.
.
.
scor(tn)
0
0
0
.
.
.
0
respectively, then
f(y) = kAy − bk
2
.
Tasks
1. Often, when defining matrices like A above, it is helpful to begin with n small enough that you can
actually look at the matrix and check that it is correct. The first cell of the M-file, titled “A small
example” demonstrates this technique. Run the cell and observe the output in the Matlab command
window. Then read the code in the cell. Notice that, since A is a sparse matrix, we are using the sparse
matrix functions speye and spdiags to create A. (Remember, if you don’t know how to use a function,
you can always use doc.)
2. Now, run the cell labeled “Part 1.” It should produce a figure like the one above. To complete this cell,
you need to define the variables A1 and b1, corresponding to A and b, in the places indicated. (Hint: you
can copy and paste from the first cell. Just make sure to rename A and b.)
Variables: A1, b1
If you look at the documentation for Matlab’s backslash operator, you will see that, when applied to an
overdetermined system, the result will be the least-squares solution. Thus, the code
% find the solution
y1 = A1\b1;
beneath your variable definitions produces the solution y that we’re looking for. If you have defined your
variables correctly, running the cell should produce the following figure.
0 0.2 0.4 0.6 0.8 1
time
-1.5
-1
-0.5
0
0.5
1
1.5
amplitude
corrupted
original
reconstructed
3
3. As you can see from the figure, the choice δ = 1 for this problem does filter some of the noise, but is still
quite rough. Try several other values for δ. Notice that δ = 1000 is too large; the smoothed signal is
essentially constant. A δ around 10 or 15 seems fairly good for this example.
4. Now, let’s try a different choice of matrix D. We’ll use the (n − 2) × n matrix
D =
1 −2 1 0 · · · 0
0 1 −2 1 · · · 0
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
0 · · · 0 1 −2 1
.
In the cell labeled “Part 2: Another small example,” define the variable D2 where indicated. Hint:
>> e = ones(5,1); full(spdiags( [e,-2*e,e], [0,1,2], 3, 5 ))
ans =
1 -2 1 0 0
0 1 -2 1 0
0 0 1 -2 1
Run the cell and make sure the output looks correct.
Variables: D2
5. In the cell labeled “Part 3,” complete the code to solve the problem (similar to “Part 1”) with the new
choice of D. Be sure the variable names for A, b and y are A3, b3 and y3 respectively. Again, try a few
choices for δ. When δ = 100, the figure should look like this:
0 0.2 0.4 0.6 0.8 1
time
-1.5
-1
-0.5
0
0.5
1
1.5
amplitude
corrupted
original
reconstructed
Variables: A3, b3, y3
Q1: Describe the y that results when δ = 105
.
