Please work through all the code and examples in this worksheet.
Things in this font are Octave commands — don't cut and paste them, type the commands yourself. Trust me, this is the better way to learn these commands.
Please do the following to show your work:
a=randn(1,100);[Note: do help randn to see the difference with the rand we have mostly been using in the past. You might do a hist(a) to see how well this conforms to what you expect after reading the help.]
Create a second random vector b with
b=randn(1,100);Plot a and b. How about plotting them on the same axes?
You can put them into a $2\times100$ matrix by [a;b], but Octave will not make a nice plot if you do plot([a;b]); instead transposing the matrix first, with plot([a;b]'), makes a much better output.
Remember the single quote ' is transpose ... in fact, it is conjugate transpose, which is slightly different in ways that mater below.
Make that nice plot of a and b.
Now, we can hide two numbers in with all the random noise in a and b by
c=3*a-5*b
How to extract the secrets 3 and -5 out of c?
Well, a and b are random, they're almost certainly nearly orthogonal, so if we take the dot product of c with a, we will mostly just recover the part of c in the direction of a. Then dividing by a${}\cdot{}$a will give the coefficient of a in the linear combination that is c.
In summary, the secrets should be (approximately) recovered by
(c*a')/(a*a') (c*b')/(b*b')since the way to do the dot product of two row vectors x and y in Octave is x*y'.
ab=[a;b]We can creeate c as above by multiplying the row vector [3 -5] by ab (areful, we're going to do a lot of left-multiplication today, because we're working with row vectors):
c=[3 -5]*aband we can compute the dot product of c with each of a and b as
c*ab'Similarly, we can compute the dot products between all possible combinations of a and b as
ab*ab'Note that this is a $2\times2$ matrix, where the top left entry is the dot product of a with itself, the bottom right entry is the dot product of b with itself, and the off-diagonal entries are the dot product of a with b. To repeat the calculations above in matrix form, we want to divide each of the entries of c*ab' by each of the diagonal entries of the matrix just computed (try entering diag(ab*ab') to be clear on what is going on here):
(c*ab')./diag(ab*ab')'This should return approximately 3 -5. Look at each part of the above expression closely to make sure you understand the parts: we first compute c*ab', which is a $1\times2$ vector; the "./" notation means element-by-element division; Octave's diag function extracts the diagonal entries of the $2\times2$ matrix ab*ab, and returns them as a column vector; the final single quote transposes this to a row vector so that the "./" operator will not return an error.
We have thus shown how to combine random vectors, and then "decode" after the fact which linear combination we used. Doing exactly these same steps below with cosines and complex exponentials will give us the cosine transform and Fourier transform respectively.
One thing worth noting here is that while this element-by-element division will approximately extract the components (3,-5), what we really want is a matrix inversion. We can invert the matrix ab*ab' via inv(ab*ab’) (instead of the -1 power...), and the following expression should return exactly 3 -5:
c*ab'*inv(ab*ab')(The ab*ab' is essentially always invertible because, as we have said in class, random matrices hardly ever singular.) This extra step will not be necessary below, because cosines and complex exponentials are so well-behaved.
plot(cos((.5:100)/100*pi*3))We can thus create a two-row array ab which we will use analogously to how we used it in the previous section, though instead of containing random rows, it will now consist of a vector of frequency 3 (as we just computed) and a vector of frequency 4, as:
ab=[cos((.5:100)/100*pi*3);cos((.5:100)/100*pi*4)];Now plot it with plot(ab'), as before.
As above, let c=[3 -5]*ab;, and perhaps plot all three vectors, with plot([ab;c]'); can you guess how to recover secret coefficients (3,-5) from this plot?
As above, (c*ab')./diag(ab*ab')'. Note that here we have recovered the coefficients exactly. Why? Because our basis vectors, the cosine functions, are or thogonal to each other: the matrix ab*ab' only has entries on its diagonal. What we have done is not quit e a Fourier transform, because we used cosine functions, instead of complex exponential functions. This is called the cosine transform.
Before we move on, it is worth noting how convenient it is that the two cosine functions we picked, cos((.5:100)/100*pi*3) and cos((.5:100)/100*pi*4) are orthogonal to each other. If we change the frequencies 3 and 4 to a different pair of numbe rs, will this remain true? Yes, if the frequencies are different integers between 0 and 99. (Where does that 99 come from? Also, remember, vectors a and b are orthogonal if their dot product is 0; due to rounding errors, you may get something like -8.9928e-015, denoting that at least the first 14 digits of their dot product is 0, and numerical error is responsible for the rest.) If you are curious, try computing the dot products ab*ab' though replacing 3 and 4 in the definition of ab above with two other numbers of your choice.
Up until now, these manipulations have not given us more than we started with: we started with a linear combination of some basis vectors, and ultimately recovered the coefficients. However, the fundamental claim of the cosine transform (and the Fourier transform) is that anyvector can be expressed in this fashion, at least, if we express it in terms of cosines of all frequencies, instead of just frequencies 3 and 4. Instead of the matrix ab with just cosines of frequencies 3 and 4, we will construct a matrix M with all the frequencies. The following is a not-very-efficient way of constructing it:
M=[]; for j=0:99, M=[M;cos((.5:100)/100*pi*j)]; end
Recall from above that we used the code (c*ab')./diag(ab*ab')' to find the linear combination of the rows of ab that produced c (provided the rows of ab were orthogonal). We compute the full cosine transform identically. If you have not changed c from above, it should still be 3 times the vector cos((.5:100)/100*pi*3) minus five times the vector cos((.5:100)/100*pi*4), and we can see this via the cosine transform, as:
plot(0:99,(c*M')./diag(M*M')')(or is stem better?). We can see (possibly after zooming in) that we have completely recovered the secret coefficients: the cosine at frequency 3 has a coefficient of 3, and the cosine at frequency 4 has a coefficient of -5.
Now that we have taken the cosine transform of something that we already know how to express via cosines, we are ready to get adventurous. Pick your favor ite 100 values and store it as f, or let f be the $1\times100$ vector such that plot(f) draws a function you like. (Try something with steps or sharp corners, it's fun.)
Plot f with plot(f) and plot its cosine transform by first computing
fcos=(f*M')./diag(M*M')';and then plotting it with plot(0:99,fcos) (or stem). Remember, the claim is that this odd pattern on screen tells you the way to combine cosines to produce the function f. This is somewhat mysterious and far from intuitive. How do smooth cosines combine to form something with corners? What is the meaning of the left half of the graph (representing coefficients of low frequency components) or the right half of the graph (representing high frequency components)? Questions like this are at the heart of the Fourier transform.
To get a sense for how these cosines can do this we should take a closer look at fcos. The claim is that fcos represents cosine coefficients, and that we can reconstruct f as fcos*M — try plotting fcos*M to convince yourself; or if you are more skeptical, compute the maximum of the absolute value of the difference between this and f with max(abs(fcos*M-f)).
In some sense this works so well that it is boring. To get a more interesting notion of what is going on here, we can break it: try using only the first 30 cosines, as in
plot(fcos(1:30)*M(1:30,1:100)The "first 30 cosines" represent the slowest-oscillating cosine functions in this expression. You could try only the fastest 30 cosines with plot(fcos(end-29:end)*M(end-29:end,1:100)). (This end is new Octave syntax. Read this page for an eplanation of this useful construction.) For most functions, the slowly-oscillating components are visually a lot more recognizable than the others. This is what many image and video compression algorithms are based on. Pl ay around with the above concepts and code until something clicks (not your mouse, in your head).
At this point it is worth taking a look at the mysterious scaling factors that we have been dividing by, diag(M*M')' — remember that we have been using such scaling factors ever since the first part of this document, when we were using random vectors. The diagonal of M*M’ stores the dot products of each of our cosines with itself; dividing by these entries, element-by-element, was just to make sure that the resulting coefficients have the right scale. Go back and review the first section if you are confused. Plot these scaling coefficients, with plot(diag(M*M')') (stem>?). It turns out that the first scaling factor is 100, and all the rest are 50. Not so mysterio us after all. In future, instead of computing this matrix product, we will use the much simpler e xpression [n n/2*ones(1,n-1)], which is equivalent but cleaner.
im=imread('filename.jpg');
Now, you have to crop it to be square, and deal with the fact that this
image matrix is 3-dimensional with the third dimension having three values
representing the red, green, and blue values in the image. You can simply
extract one of the colors, or combine them, or something like
im=sum(im(1:500,end-500:end,:),3);— tinker around a bit!
Now that you are dealing with only one value of image intensity, not three of colors, to picture this you have to tell Octave to use a "colormap" in which 0 is black and large numbers are white:
colormap('gray')
imagesc(im)
axis image
(If you like, also do colorbar to get a visualization of which colors
correspond to which numbers.)
We can now try some cosine transforms on this image. The image here was $500\times500$ (yours will likely not be! adjust as necessary in the following), which means we should recalculate our transform matrix M:
n=500;M=[]; for j=0:n-1, M=[M;cos((.5:n)/n*pi*j)]; endAs mentioned at the end of the last section, instead of computing M*M' each time we compute coefficients, we can instead just use the vector [n n/2*ones(1,n-1)] — assuming you do not change the variable n in the interim. Thus to take the cosine transform of all the rows in the image, we multiply the image array by M', and then divide each row element-wise by [n n/2*ones(1,n-1)] (recall that multiplying a matrix A by a matrix B can be thought of as multiplying each row of A by the matrix B). In order to rescale the elements of each row, we could use a for loop, but it might be conceptually easier to just replicate the scaling factor [n n/2*ones(1,n-1)] so that there is one copy of it for each row. One clever way to do this is to multiply ones(n,1) by it. Thus the cosine transform of the rows of im is
imR=(im*M')./(ones(n,1)*[n n/2*ones(1,n-1)]);Feel free to look at our result now, with imagesc(imR), but we are not done yet: we also need to cosine-transform each column. The easiest way to do this is to flip the rows and columns of imR, do the cosine transform on each new row, and then flip rows and columns back:
imcos=((imR'*M')./(ones(n,1)*[n n/2*ones(1,n-1)]))';You can look at this with imagesc(imcos), though you may need to zoom in on the top left corner before you see anything: only the low-frequency cosine comp onents of an image have much action. To see something about the whole structure of the cosine transform, you might have to squash the range, for example, by displaying the fourth-root of the absolute values:
imagesc(abs(imcos).^.25)We can recover the original image by multiplying by M, then transposing and multiplying by M again (to transform the columns), and then transposing again:
imagesc(((imcos*M)'*M)')However, it might be interesting to modify the cosine transform before we invert it. We will use imcosm to play around with modifications. We can, for example, take the first $30\times30$ cosine components and make them 0:
imcosm=imcos; imcosm(1:30,1:30)=0; imagesc(((imcosm*M)'*M)')This shows us only the high-frequency, rapidly-oscillating parts of the image. Alternatively, we could see only the components that are oscillating slowly in their first coordinate (vertically, in Octave), with:
imcosm=imcos; imcosm(31:end,:)=0; imagesc(((imcosm*M)'*M)')Play around with this for a bit. For example, one other odd effect is if we negate the small cosine components:
imcosm=imcos; imcosm(1:30,1:30)=-imcosm(1:30,1:30); imagesc(((imcosm*M)'*M)')One bit of terminology, so that you can speak with the rest of the world about what you have learned: the process of expressing a signal in terms of cosines as we have been doing is called the inverse (discrete) cosine transform, which for our vector c we computed as (c*M')./diag(M*M'), while the process of taking coefficients and computing the corresponding linear combination of cosines is called the (discrete) cosine transform, which we computed by just multiplying by M.
Now that we have seen the cosine transform, the Fourier transform can be seen as just a minor variation. The Fourier transform uses complex numbers, however, which can add some headaches. In particular, one thing to be careful of is that the single-quote operator we have been using to transpose matrices in Octave actually computes the complex conjugate of the transpose, which turns out to be fairly useful. But if we want just the transpose, add a dot in front of the quote, as in M'.
While the cosine transform involved expressing a signal as a linear combination of cosines, $\cos(kx)$, for different $k$, the Fourier transform involves expressing a signal as a combination of the complex analog of this, $e^{-ixk}$ = $\cos(xk) - i\sin(xk)$ for different $k$. For example, the function that oscillates 3 times around the unit circle between 0 and 100 is:
z=exp(-1i*(0:99)/100*2*pi*3)These expressions, where we change 3 to any value between 0 and 99, are known as the Fourier basis.
We can now plot the real and imaginary parts of $z$ with plot(0:99,[real(z);imag(z)]') (real is in blue, imaginary is in green), or we can plot them in the complex plane even more simply, with plot(z), which is equivalent to plot(real(z),imag(z)). Since our function goes around the unit circle three times, if you zoom in you should be able to see that the circle is actually drawn 3 times.
It is worth noting that the negative sign inside the exponent is a convention, and different definitions of the Fourier transform will involve different signs, possibly multiplying the answer by $n$ or $\sqrt{n}$ or something involving $\pi$. None of these variants fundamentally change what is going on.
In analogy with how we constructed the matrix M for the cosine transform above, we will now construct the corresponding matrix for the Fourier transform (in general, this document keeps reusing code with only minor variations; try to keep track of the patterns, and the intuition behind these patterns):
n=100;F=[]; for j=0:n-1, F=[F;exp(-1i*(0:n-1)/n*2*pi*j)]; endThe functions in the rows of this matrix are the Fourier basis. Comparing this with our construction of M above, we see there are a few minor differences: the range inside the exponential starts at 0 instead of .5, and there is an extra factor of 2, but other than that, the Fourier transform is very similar to the cosine transform.
To do a Fourier transform, that is, given 100 coefficients, computing the corresponding linear combination of the basis functions, can be done as above simply by multiplying by F . To convince ourselves that we are doing the right thing, we can create a random vector r, and compare our Fourier transform, to Octave's built-in fft function:
r=rand(1,100); rf=r*F; max(abs(fft(r)-rf))Computing the inverse of this, namely, figuring out how to express a function as a linear combination of the Fourier basis, is the same as above: given rf , we compute:
r2=(rf*F')./diag(F*F')';This should now be very close to r, which we can verify with max(abs(r-r2)). Actually, if you look at the scaling factors diag(F*F')' you will see that they are all 100, namely n, so we can simplify our expression for the inverse Fourier transform to just: multiply by F' and divide by n.
While in the previous section we took the cosine transform of our image, we can now take the Fourier transform. First construct the transform matrix for n = 500:
n=500;F=[]; for j=0:n-1, F=[F;exp(-1i*(0:n-1)/n*2*pi*j)]; endAnd then transform the rows, transpose it (with .' now since the values are complex), transform again, and transpose back:
imF=((im*F).'*F).';We can display the results, as above, with imagesc(abs(imF).^.25) where the .25 power is to compress the range to where we can see it. Note that, unlike for the cosine transform, there is now action in all four corners of the image. This is because, while the top left corner of the image still corresponds to low frequencies in the $x$ and $y$ directions, the bottom right corner corresponds to low negative frequencies, and the top right corner corresponds to low negative frequencies in the $x$ direction and low positive frequencies in the $y$ direction, etc. It makes intuitive sense to "circularly shift" the image so that all four corners are plotted together in the center of the image, which we can do with imagesc(circshift(abs(imF).^.25,[199,199])) where the parameter [199,199] means to shift by 199 (which is $n/2$) in each dimension. To get a nice-looking figure, change the colormap to colormap('hsv').
n=4;F=[]; for j=0:n-1, F=[F;exp(-1i*(0:n-1)/n*2*pi*j)]; end; FIn row $j$, column $k$ (indexing from 0), the value will be $e^{-2\pi ijk/n}$, which, if we let $\omega$ be defined to be $e^{-2\pi i/n}$, means that the $(j, k)$th entry of F is just $\omega^{jk}$. The exponents form a multiplication table! Further, the matrix F is symmetric — swapping $j$ and $k$ (swapping rows and columns) does not change anything. The Fourier transform is multiplication by a matrix consisting of $\omega$ raised to powers that form a multiplication table. Remember that $\omega$ itself is what is called a primitive $n$th root of unity, meaning that if you raise it to the $n$th power, it will be 1. (Primitive means means that no power smaller than $n$ will take $\omega$ to 1.)
Since multiplication by F gives the Fourier transform, and multiplication by F'/n gives the inverse Fourier transform, it turns out that the Fourier transform and its inverse are really almost the same thing: the only difference between F and F' is that F' is the complex conjugate of F (recall that for complex matrices, the single quote operator in Octave takes the complex conjugate as well as transposing the matrix; Since F is symmetric, however, transposing it does not change it, and the single-quote operator's only effect is to take the conjugate of F).
Thus in Octave, instead of using the inverse Fourier transform function ifft, we could instead just take the complex conjugate of the fft of the complex conjugate of our vector, and divide by the length of the vector. We can compare these two approaches for a random complex vector r:
n=100; r=rand(1,n)+1i*rand(1,n); max(abs(ifft(r)-conj(fft(conj(r)))/n))
conv([1 1 1 1],[1 1 1])The basic fact underlying this section (and in fact a lot of the practical usefulness of convolution) is that convolution becomes element-by-element multiplication in Fourier space.
With convolution, you can take an image and make two copies of it, separated horizontally by 100 pixels, by convolving the image with the vector [1 zeros(1,99) 1] — try it. We can do this in Fourier space as follows (remember to reset colormap('gray') if you had changed the colormap earlier). First we will create the convolution kernel, which has two 1's in it, separated by 100 horizontally, and has the same size as the image:
kern=zeros(500,500); kern(1,1)=1; kern(1,101)=1;Then we can convolve this with im and plot the results as:
imagesc(ifft2(fft2(im).*fft2(kern)))(The "2" here tells Octave to do FFT along both dimensions!) You should now be looking at two copies of your image, 100 pixels horizontally shifted from each other. One slightly strange thing is that the right edge of the image wrapped around to the left, for the copy that we shifted right by 100 pixels — we might have hoped this operation would expand the image to be 100 pixels wider in the horizontal direction, but of course the Fourier transform does not change the size of the arrays it operates on. (On
Try coming up with your own convolution kernels:
kern=zeros(500,500); kern(1:10,1:10)=1;imagesc(ifft2(fft2(im).*fft2(kern)))The above kernel blurs the image by taking averages of 10-by-10 blocks of pixels (though note how the bottom of the image wraps around and bleeds into the top).
kern=zeros(500,500); for j=1:10,kern(j,j)=1;end;imagesc(ifft2(fft2(im).*fft2(kern)))The above kernel blurs the image diagonally.
kern=zeros(500,500); kern(1,1)=1; kern(1,2)=-1;imagesc(ifft2(fft2(im).*fft2(kern)))Negative numbers lead to more interesting effects: the above kernel computes the difference between horizontally adjacent pixels, effectively becoming a "vertical edge detector".
kern=zeros(500,500); kern(1,1)=1; kern(2,1)=-1;imagesc(ifft2(fft2(im).*fft2(kern)))Switching things slightly leads to a vertical edge detector. Or, for something a bit stranger, if you looking at a particular feature $\Ff$ in your image, we can build a "$\Ff$-detector" by using convolution to try to match a part of $\Ff$:
kern=zeros(500,500); feat=im(199:-1:100,379:-1:320);kern(1:100,1:60)=feat-mean(feat(:)); imagesc(ifft2(fft2(im).*fft2(kern)))One of the most intriguing features of convolution via Fourier transforms is that this methodology also lets you invert the operation. Suppose you know an image has been blurred in a certain way (perhaps motion-blurred in a certain direction?) then, instead of multiplying by the Fourier transform of the kernel, we can divide by it to invert the operation. In some sense, you can invert blur filters to get sharpening filters. Be careful though, if any of the Fourier coefficients of your kernel are 0 or small, then dividing by it could lead to divide-by-zero errors, or numerical blowup in other ways.
One particularly useful application of this is to solve differential equations: derivatives are in some sense convolutions (we computed "cheap" derivatives above to detect horizontal and vertical lines), so if $f^\prime$ and $f^{\prime\prime}$ can be expressed as convolutions of certain kernels with $f$ , (and of course $f$ itself is the convolution of the trivial kernel 1 with $f$ ), then any linear differential equation can be represented via convolution as $conv(f, K) = g$, and we can solve the differential equation by taking the inverse-convolution of $g$ with $K$ to find $f$, via some variant of $f = {}$ifft$($fft$(g)/$fft$(K))$.
Recall the matrix F that we had constructed, for various sizes, to help us compute the Fourier transform. We need to understand F well enough so that we can multiply by F in a divide-and-conquer fashion. n) time. We will construct F in an easier way than we did above — now that we understand what F means, we can cheat and use Octave's built-in fft function:
F=fft(eye(8))The fft function transforms each of the rows of this matrix. Essentially, the first row of F tells us how Octave transforms entries in the first column of a vector; the second row of F tells us how Matlab transforms entries in the second column of a vector, etc., which is exactly how we defined F above. (As a side comment, the fft function actually transformed each column of its input, not each row, but since F is symmetric, this is the same thing.)
The first thing we note about F is that the submatrix consisting of rows 1,3,5,7 (indexed from 1) and the first four columns, is the same as the submatrix of F consisting of the last four columns of these rows:
F(1:2:end,1:4) F(1:2:end,5:end)We can easily prove this fact by recalling the definition of F from the "Structure of the Fourier transform matrix" section, namely that for $\omega = e^{-2\pi i/n}$ being an $n$th root of unity, the $(j, k)$th entry of F, indexing from 0, is $\omega^{jk}$. The claim is that for even $j$, and $k < \frac{n}{2}$, $\omega^{jk}=\omega^{j(k+\frac{n}{2})}$. Since $j$ is even, $j\frac{n}{2}$ in the exponent is a multiple of $n$, and since $\omega^n=1$, this term will vanish, leaving just $\omega^{jk} as claimed.
This is progress towards our goal of a divide-and-conquer algorithm for multiplying by F, that two parts of the multiplication by F are actually the same. But what we would really like is for them to be related to a Fourier transform of a smaller size. As it turns out, we are in luck: the matrix F for Fourier transforms of size 4 is
fft(eye(4))And this matrix is identical to the two sub-matrices of the $8\times8$ transform matrix we computed previously. Thus we have divided the $8\times8$ Fourier transform matrix into four $4\times4$ submatrices (though we were a bit tricky with the rows), and we have looked at two of these submatrices and found that they are the $4\times4$ Fourier transform matrix, which means our divide-and-conquer algorithm can solve all of these recursively.
Now we need to worry about the two other portions of F, the even rows (or odd rows, if you are thinking in terms of 0-indexing):
F(2:2:end,1:4) F(2:2:end,5:end)There are two disasters here: first, these two matrices are not equal, and second, neither of them equals the $4\times4$ Fourier transform matrix.
The first disaster, though, is easily averted: the second matrix is just the negation of the first matrix. And the second disaster turns out to be as benign: each row of these matrices is a fixed multiple of the corresponding row of the $4\times4$ Fourier transform matrix, meaning that to compute the product of a vector with one of these matrices, we just multiply the vector by the $4\times4$ Fourier transform matrix, and then multiply each entry in the answer by the appropriate scale. Perhaps you can figure out by now why the $n/2 = 4$ scaling coefficients should be $(1, \omega, \omega^2 , \omega^3)$ for $n = 8$, with the natural generalization for other $n$.
We will write a function myfft to implement the full divide-and-conquer Fourier transform algorithm (for vectors whose length is a power of 2). Save the following code as myfft.m
function x=myfft(x)
n=length(x);
omega=exp(-1i*2*pi/n);
scale=omega.^(0:n/2-1);
if n>1,
even=myfft(x(1:2:end));
odd=myfft(x(2:2:end));
x=[even+scale.*odd, even-scale.*odd];
end
endfunction
Spend some time and try to figure out how this code corresponds to the above
explanation; or write your own code if that would be easier, or check
Wikipedia. We can check that the code is doing the right thing by comparing
it with the built-in fft, on a vector of length, say, 1024:
r=rand(1,1024); max(abs(fft(r)-myfft(r)))