Math551 - lab10 -Solved

Goal: In this lab, you will learn to use Matlab’s backslash operator for signal denoising.

Getting started: You will need to download the file m551lab10.m and put it in your working directory.


Variables: x1, A, q
At times you will also be asked to answer questions in a comment in your M-file. You must format your text answers correctly. Questions that you must answer are shown in gray boxes in the lab. For example, if you see the instruction

Q7: What does the Matlab command lu do? your file must contain a line similar to the following somewhere

% Q7: It computes the LU decomposition of a matrix.

The important points about the formatting are

•    The answer is on a single line.

•    The first characters on the line is the comment character %.

•    The answer is labeled in the form Qn:, where n stands for the question number from the gray box.

If you do not format your answers correctly, the grading software will not find them and you will lose points.

Basic Ideas
Consider the following figure.

 

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

 s(t1)

s(t2) x =  ...            and         xcor.

s(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 Rn, define

n−1

f(y) = ky − xcork2 + δ2 X(yi+1 − yi)2.

i=1

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

−1

0



D =  ...





0

Then
1

−1

...

···
0

1

...

0
···

···

...

−1
0

0

...

1
so that
 y2 − y1   3 − y2  y

Dy =         ...   .



yn − yn−1
f(y) = ky − xcork2 + δ2kDyk2.
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 the (2n − 1)-vector b as

                                     A        and     b  

respectively, then

f(y) = kAy − bk2.

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.

 

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

                                                                                         1     −2       1       0       ···      0

                                                                                         0      1       −2      1       ···      0

                                                                             D = ...    ...        ...       ...       ...    ....

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:

 

Variables: A3, b3, y3
Q1: Describe the y that results when δ = 105.