DFT part 2: randomising data and filtering

I. Cnop
Vrije Universiteit Brussel
icnop@vub.ac.be

Prerequisite: Building the DFT

See the appropriate notebook

Sum of a randomised progression.

Introduction

In the previous notebook DFTBuildup we built a wheel with spokes at multiples of  2π/n, called directions.  Since multiplying angles corresponds to powers of points, adding all directions forms a simple geometric series that adds op to zero. The proof is simple: adding one by one all line segments (starting with the spoke at the vertex 1) gives a regular polygon that closes at the origin. If we hang the spoked weel from a thread it is in perfect balance.

We will now modify the lengths of the spokes slightly and see what happens.

Recall the shorthand for directions in the plane

directions[
        n_] := (E^(I Range[n] 2 π/n ))

coords[z_] := {Re[z], Im[z]}

Add some random value to lengths

coeffList = Table[1 +   .2 Random[ ], {13}]

n =  Length[coeffList]

Multiply

Map[coords, coeffList directions[n] // N];

ListPlot[%, AspectRatio -> Automatic, 
      PlotRange -> {{-2, 2}, {-2, 2}}]

Add

Making partial sums by adding these points is done wit a built-in  Mathematica  command

points = FoldList[Plus, {0, 0}, %%];

ListPlot[points, PlotJoined -> True, 
      AspectRatio -> Automatic, 
      PlotRange -> {{-4, 5}, {-2, 5}}]

The final result

We finally get back close to the origin, but not quite.

Experiment

Repeat this for other choices of random modifications of the constant list above . Here is one that averages near one:

coeffList = Table[ .8 +   .4 Random[ ], {13}]

n =  Length[coeffList]

Map[coords, coeffList directions[n] // N];

ListPlot[%, AspectRatio -> Automatic, 
      PlotRange -> {{-2, 2}, {-2, 2}}]

points = FoldList[Plus, {0, 0}, %%];

ListPlot[points, PlotJoined -> True, 
      AspectRatio -> Automatic, 
      PlotRange -> {{-4, 5}, {-2, 5}}]

Maybe it is even more apparent if we add segments

Map[Line, 
        Transpose[{Table[{0, 0}, {n}], 
            Map[coords, coeffList directions[n]]}]];

Show[Graphics[%], AspectRatio -> Automatic, 
      PlotRange -> {{-4, 5}, {-2, 5}}]

The endpoints of the segments constitute the vertices of a slightly irregular  n -gon. 
Joining all segments as matchsticks one after the other draws the bigger   n - gon in the previous Graphics cell.

The reason that we end up close to the origin is that errors going in opposite directions level out. 
This is quite similar to the rolling of dice (where sums tend to a multiple of 3.5) or tossing a coin if we replace heads by 1, tails by -1 and repeating the toss.

It can be explained by the analogue of a drunk person :by stumbling over some random distance (with distances averaging to 1) in each of the n directions   2πj/n    he will not get far.

List programming for a sum

Obtaining the final sum is easiest with a   Dot  product command: try

afew = {a, b, c, d}
m = Length[afew]
counts = Range[m]
afew counts
afew . counts

We have performed

coeffList = Table[ .8 +   .4 Random[ ], {13}];

n =  Length[coeffList];

coeffList . directions[n]

The fact that these values are close to the origin is enhance if we divide by  the square root of  n  , which we did to obtain the discrete Fourier transform.

Discrete Fourier transform of a noisy list.

Vectors with arbitrary lengths

We recommend that this buildup be first performed with a list of positive numbers for lengths. 
However it works as well with negative lengths; in this case the spokes representing vectors are sticking backward from the direction.
Her is an irregularly oscillating list

coeffList = 
        Table[(((-1))^j) ( .8 +   .4 Random[ ]), {j, 
            14}];

ListPlot[coeffList]

coeffList // TableForm

For later reference, here is their sum:

Plus @@ coeffList

n =  Length[coeffList]

Map[coords, coeffList directions[n] // N];

We now plot the endpoints of the vectors with lengths in  coeffList

Remark: since some numbers are negative, some of these endpoints lie on the other side of the origin, back from the original ray. 
See this phenomenon in the following graphic:

ListPlot[%, AspectRatio -> Automatic, 
      PlotRange -> {{-2, 2}, {-2, 2}}]

Adding vectors one by one

produces a kind of walk from the origin:

points = FoldList[Plus, {0, 0}, %%];

ListPlot[%, PlotJoined -> True, AspectRatio -> Automatic, 
      PlotRange -> {{-2, 2}, {-2, 2}}]

But only half of the spokes are shown in the following

Map[Line, 
        Transpose[{Table[{0, 0}, {n}], 
            Map[coords, coeffList directions[n]]}]];

Show[Graphics[%], AspectRatio -> Automatic, 
      PlotRange -> {{-2, 2}, {-2, 2}}]

We can add arrows using the

<< Graphics`Arrow`

package, to see overlapping spokes.

Selection of the k-th harmonic

Start by counting angles by  multiples   of   a fixed number    k   :

mixDirections[n_, 
        k_] := (E^(I k Range[n] 2 π/n ))

We are in fact using multiplication in  (Subscript[U, n])

Table[Mod[k Range[13], 13], {k, 13}] // TableForm

We will use this way of counting to make new sums.

Selection of the 2-nd harmonic

Map[coords, coeffList mixDirections[n, 2] // N];

ListPlot[%, AspectRatio -> Automatic, 
      PlotRange -> {{-2, 2}, {-2, 2}}]

makes two turns around the origin and is grouping even-numbered and odd-numbered lengths from coeffList before adding:

points = FoldList[Plus, {0, 0}, %%];

ListPlot[%, PlotJoined -> True, AspectRatio -> Automatic, 
      PlotRange -> {{-2, 2}, {-2, 2}}]

Both end result and partial sums are strikingly different. Again we may also add lines

Map[Line, 
        Transpose[{Table[{0, 0}, {n}], 
            Map[coords, coeffList mixDirections[n, 2]]}]];

Show[Graphics[%], AspectRatio -> Automatic, 
      PlotRange -> {{-2, 2}, {-2, 2}}]

Selection of the 3-rd harmonic

Counting by three still gives a new result:

Map[coords, coeffList mixDirections[n, 3] // N];

ListPlot[%, AspectRatio -> Automatic, 
      PlotRange -> {{-2, 2}, {-2, 2}}]

points = FoldList[Plus, {0, 0}, %%];

ListPlot[%, PlotJoined -> True, AspectRatio -> Automatic, 
      PlotRange -> {{-2, 2}, {-2, 2}}]

Map[Line, 
        Transpose[{Table[{0, 0}, {n}], 
            Map[coords, coeffList mixDirections[n, 3]]}]];

Show[Graphics[%], AspectRatio -> Automatic, 
      PlotRange -> {{-2, 2}, {-2, 2}}]

Selection of the 4-th harmonic

Map[coords, coeffList mixDirections[n, 2] // N];

ListPlot[%, AspectRatio -> Automatic, 
      PlotRange -> {{-2, 2}, {-2, 2}}]

points = FoldList[Plus, {0, 0}, %%];

ListPlot[%, PlotJoined -> True, AspectRatio -> Automatic, 
      PlotRange -> {{-2, 2}, {-2, 2}}]

Map[Line, 
        Transpose[{Table[{0, 0}, {n}], 
            Map[coords, coeffList mixDirections[n, 4]]}]];

Show[Graphics[%], AspectRatio -> Automatic, 
      PlotRange -> {{-2, 2}, {-2, 2}}]

The 0-th harmonic

Is just the sum of all entries in coeffList

These harmonics are all close to the harmonics for the unrandomised list since the contribution of the random part is small (an could even be neglected).

Fourier[coeffList] // TableForm

Fourier[Table[(((-1))^j), {j, 14}]]

error = % - %%

Map[coords, error];

ListPlot[%, AspectRatio -> Automatic, 
      PlotRange -> {{-2, 2}, {-2, 2}}]

Animations are provided for the rearrangements of spokes

Do[
Show[
        Graphics[
Map[Line, 
            Transpose[{Table[{0, 0}, {n}], 
                Map[coords, 
                  coeffList mixDirections[n, k]]}]]
], 
        AspectRatio -> Automatic, 
        PlotRange -> {{-2, 2}, {-2, 2}}], {k, n}]

Run these slowly.

Filtering to obtain the original coefficient List.

We want to eliminate the effects of the randomness before computing th original list bij applying the inverse Fourier transform. To establish the kind of uncertainty we want to eliminate from the Fourier transform we have to know some standard profiles for these.

The original list

coeffList = Table[(((-1))^j), {j, 14}];

ListPlot[coeffList]

coeffList // TableForm

For later reference, here is their sum:

Plus @@ coeffList

n =  Length[coeffList]

Map[coords, coeffList directions[n]];

ListPlot[%, AspectRatio -> Automatic]

We see seven points , each  counted twice.

Fourier[coeffList];

ListPlot[Map[coords, %], AspectRatio -> 1, 
      PlotRange -> All]

We see one dot to the left of the diagram and and all other harmonics are zero.

The modified list

coeffList = 
        Table[(((-1))^(j )) ( .9 +  .2 Random[ ]), {j, 14}];

ListPlot[coeffList]

coeffList // TableForm

For later reference, here is their sum:

Plus @@ coeffList

n =  Length[coeffList]

Map[coords, coeffList directions[n]];

ListPlot[%, AspectRatio -> Automatic]

We see fourteen points that do not overlap any more.

Fourier[coeffList];

ListPlot[Map[coords, %], AspectRatio -> 1, 
      PlotRange -> All]

We see 1 value far out, and most other values close to the origin (in fact within a distance  ABS  less than  .1)

We will eliminate the latter before inverting Fourier:

clean[list_] := 
      Table[If[Abs[list[[j]]] <  .2, 0, list[[j]]], {j, 
          Length[list]}]

InverseFourier[clean[Fourier[coeffList]]]

ListPlot[%]

This is a nice alternating (constant) sequence with values near 1 and -1.

Conclusion

By eliminating unexpected inaccuracies in the list of Fourier of a randomised list we can eliminate the randomness from the original list.
It takes some experience to recognise the unwanted inaccuracies and this expertise can be obtained by computing several standard  DFT  images.

The corresponding  Mathematica  commands are

Fourier

InverseFourier

Experiment

Try this forthe following coeffList and its modification:

coeffList=Table[(((-1))^j) +  .5 (((-1))^Quotient)[j, 4], {j, 16}];
ListPlot[
      coeffList]

n =  Length[coeffList]

Map[coords, coeffList directions[n]];

ListPlot[%, AspectRatio -> Automatic]

Fourier[coeffList];

ListPlot[Map[coords, %], AspectRatio -> 1, 
      PlotRange -> All]

coeffList=Table[(((-1))^j) ( .8 +  .4 Random[ ]) +  .5 (((-1))^Quotient)[j, 4], {j, 16}];
ListPlot[coeffList]

n =  Length[coeffList]

Map[coords, coeffList directions[n]];

ListPlot[%, AspectRatio -> Automatic]

Fourier[coeffList];

ListPlot[Map[coords, %], AspectRatio -> 1, 
      PlotRange -> All]

What cleaning do you expect here? 
Try to remove everything at a distance smaller than, say .1, and do the inverse Fourier transform.

clean[list_] := 
      Table[If[Abs[list[[j]]] <  .1, 0, list[[j]]], {j, 
          Length[list]}]

InverseFourier[clean[Fourier[coeffList]]]

ListPlot[%]

The threshold value can be reset without much effect: we can clean with a coarser brush without affecting the final result:

clean[list_] := 
      Table[If[Abs[list[[j]]] <  .4, 0, list[[j]]], {j, 
          Length[list]}]

InverseFourier[clean[Fourier[coeffList]]]

ListPlot[%]

We do not recover the unmodified coeffList, but there is more regularity in these numbers. 
Left and right pairs ( pairs 1, 3 and 9, 11, pairs  4, 6 and 12,14  ,...) now align in a similar way. 
Compare this with the original input where this is not the case:

Large numbers of data give better results.

The same reasoning can be performed when  n  gets larger. 
A law in probability states that the larger the number of data, the smaller the effect of a random error on these data.

Task: filtering noise from sound.

We are going to repeat the construction for a very noisy (highly randomised) input of about 1000 numbers.
The set of values is wavelike, and corresponds to sampling of a basic sound. Sounds may be considered as combinations of sine functions, and a computer representation (think of MP3) is such a sampling. We are going to see what happens if we corrupt the original sound by random noise. In the example this noise is very heavy, with an intensity as loud as the original sound).

n = 1000.

coeffList = 
        Table[ Sin[2 j Pi/n] + 1.  Random[ ], {j, 1000}];

ListPlot[%]

Fourier[coeffList];

ListPlot[Map[coords, %], AspectRatio ->  .05, 
      PlotRange -> All]

This yields a cloud around the origin of small values, together with a couple of dots far from the origin.

Fourier[Table[ Sin[j Pi/n], {j, 1000}]];

ListPlot[Map[coords, %], AspectRatio ->  .05, 
      PlotRange -> All]

The latter has a tighter concentration around the origin. 
This leads us to believe that we could replace the cloud in the first graph by the origin to remove most of the random noise.
We can experiment to see with which value of "small" this happens

small =  .65

error = 
        Fourier[coeffList] - 
          Fourier[Table[ Sin[j Pi/n], {j, 1000}]];

Map[coords, error];

ListPlot[%, AspectRatio -> Automatic, 
      PlotRange -> {{-2, 2}, {-2, 2}}]

filter = 
        Table[If[Abs[Fourier[coeffList][[j]]] < small, 
            0, Fourier[coeffList][[j]]], {j, 
            Length[coeffList]}];

InverseFourier[filter];

ListPlot[%, PlotRange -> {-1, 2}]

And compare this with the original, adding the average of all random errors

Table[ Sin[2  j Pi/n] +  .5, {j, 1000}];

ListPlot[%, PlotRange -> {-1, 2}]

Experiment with different values of "small" to see when the noise dissapears.

small =  .5

small =  .6

small =  .7

small =  .75

Conclusion

The DFT and the above filtering technique is the basis of many applications in image manipulation and signal processing.
    Image enhancing and medical imaging use DFT.
    Scanners in medicine,radar images,speech technology,image compression and CD records would not exist without this technique.