What is a Discrete Fourier Transform? | Week 14 | MIT 18.S191 Fall 2020 | Grant Sanderson

For this last

week

I wanted to give a sort of one-off lecture on

what

might seem like a completely different topic, which is the fast Fourier

transform

, but Fourier

transform

s are one of those topics that really permeates all of mathematics and, as a result, It also permeates a lot of computing if you look back at the first unit we did, we were talking about image processing and one of the topics there was convolutions. Now in that lecture I talked very briefly about the computational complexity of convolutions and the fact that if you use Fourier transforms in the right way, you can make the computational complexity much smaller than you might think and then in the epidemic unit, we were talking about the type of data that is best represented with a time series.

Fourier transforms are about taking time series data, but instead of representing it in terms of frequencies in a way that I will show you in a moment and this can often make calculations faster or sometimes more intuitive and simply gives them a different vision. language with which to analyze that type of data and then later in the climate module we talk about all kinds of different differential equations to include the diffusion equation and Fourier transforms and in particular something called Fourier series can be very helpful when it comes to actually solving Those analytically aren't something we've talked about, but hopefully they give you an idea that this really permeates all kinds of different things, even if they seem disparate.

Now there are many different variations on the basic theme of a Fourier transform, so you will hear terms. like the fast 48 transform, the

discrete

Fourier transform, the Fourier series, things like that, so the word Fourier transform generally refers to something that is done with continuous signals in a Platonic ideal where you have infinitely many points, of course , in computers we only have a finite number of mini, effectively, a function is It is treated as a list of values, so we apply

what

is called a

discrete

transform or, sorry, a discrete Fourier transform and that is in the one that you and I are going to enter first today.

I'm just going to show you what it does. and then we'll talk about what the fast Fourier transform is, so I think the best context to describe Fourier transforms is with sound and musical notes, so I recorded myself humming a couple of different notes and then placing them on top of each other . So I have four different tracks here and all together this is what they sound like and they're basically broken down into individual sounds each of which is a note so for example there's the bottom note and then I saved it as its own file and then similarly. there's the note which is a little bit of that and I save it as its own file and if we go directly to our notebook here in the same directory where this notebook is located, I've saved five different sound files, so one of them plays alone . contains um, the track for the lowest note and then similarly, it sounds two, it sounds three, it sounds four, each one has a track that is just one of those notes and then it sounds one, two, three, four, contains the track with all of them placed on top of each other, so the first thing we need to do is read this data to analyze it, so I already imported and added all the different packages that we are going to need for this laptop and one of them is wav or such Maybe I should say wave.

I never actually knew if you're supposed to pronounce it dot wave or dot wav, but in any case, if we take wav read and call it in one of these files, I'll just pull out the first file. it gives us a lot of data we can work with but it looks like a mess and if we take a look at the documentation for web reading it will tell us that it knows the format of the inputs we want and if we scroll down it tells us what returns and returns a four tuple where the first term is an array containing the waveform samples, so the samples are actually the heart of the sound data which is what we need to analyze and then the rest of in It's actually metadata about the file, so what I would really care about here is that let's say we give it a name, we could call it, you know, sound and if I take out just the first element of that sound, then it's an array and it looks like. like it's a pretty long array, it's got 90,000 different elements, so if I plot this, we should get a graph that looks like a waveform, and since it's just this constant volume hum, what we're seeing is just the fact that it fades out, it stays at a constant volume and then it fades out and if we zoomed in on some of it to make it look like it actually starts to get to around uh 20,000 as an input, so we could set our x limits to You know, 20 000 and it goes up to 22,000, let's say and see how it looks.

Oh, it says x lim is not defined correctly because I'm setting it as an argument, um, or it's x limbs, okay, yeah, so x lens so what. What we see is that it looks like a wave, it's not a perfect sine wave because you know a human voice is not a synthesizer, but it produces some type of wave with a certain frequency and if we looked at a different one of the sounds, um. We may see a different frequency, so let's go ahead and prepare to do that. Here I'm taking the sound as one of those files, but let me stream it to all the files and call it something like sounds so we have several different ones and then what I'll do is just take out the signal that we really care about so that each of the signals be, you know, the first element of the sound for sound in sounds, my list of sounds and then here the thing.

What I'm plotting is one of the signals, so maybe we'll say the signals are great, but now I've removed all the different signals, so I can also do it for the second sound or the third sound, which were those hums at different frequencies. . so if I plot the second signal and it has the same x limits we see that it looks very similar but has a higher frequency and similarly the third signal is an even higher frequency it oscillates more times per second and the fourth signal is an even higher one and then the last one, which are all combined on top of each other, it seems more complicated, it's not as simple as a single sine wave.

Now before we call the Fourier transform, there is a slight awkwardness where these signals meet. read in um as two different arrays that have to do with the file tracking a left channel and a right channel, but really for us it will be better if we only have one dimensional data to deal with for each of these signals. that I'm taking out, I'm going to specify that we only want one column, so I'm going to index the entirety of a column, but saying that we only want the first one, so if we do that, when we graph You know, we just see that there are a series of data you're trying to plot for us, which will make things a little easier.

Now that we have all our signs, we have an idea of ​​what they look like. We are going to call this function fft and if we look at its documentation, it refers to the fact that a one-dimensional fft calculates the one-dimensional discrete Fourier transform and defines it with a certain formula that we are going to delve into with an animation in a moment, so the phrase discrete Fourier transform refers to the actual mathematical operation we are doing and the term fast Fourier transform refers to a specific algorithm for implementing the discrete Fourier transform that makes it actually very fast, faster than you think. would be expected.

We'll talk about it a little bit later, but now just to show what it does if we call this function on one of the signals, what we get is an array of complex values ​​and it will actually be the same size as the original signal so yeah we call this you know ff so we'll call it fft signal right and then I ask what is the length of this length of my fft signal ah give me just a moment um it's about 90,000 and that's actually about the same as the length of the signal itself so if we were calling signal length or I guess I just have a series of signals so I index them to be the same length and since I'm dealing with um you already know the list of different signals let's move on go ahead and let's say we want to stream this fast

fourier

transform across all the different signals we have so we can see what it looks like for the five different sound files, so the ffts signals are going to be this streaming type of thing. and if I, if I just trace one of them to take the signal, maybe I'll draw the first one and I'm not going to try to trace everything.

I'm only going to plot the first thousand of them because before when I tried to plot everything it really slowed everything down and the reason is because by default it tries to show you something in the complex plane because the list of values ​​that we're getting are complex numbers, so which maybe makes it harder to think about them. Now, in a moment, we'll talk about what these outputs are on the complex plane and why complex numbers have come into play and how to interpret them, but for now, I just want to show you what this does to our sound files.

I will say that the most salient information for each of those complex number outputs that we're going to look at is simply its size, its total magnitude, how far it is from the origin. So if we call absolute value and we broadcast this to all of these different complex numbers and that's what we plot, what we end up getting is this very nice graph where we have a bunch of different peaks at different heights, so what does each one represent? ? of these peaks corresponds to a frequency actually, so, for example, the first sound file that outputs something that's almost like a pure sine wave, and if we go back, it's the one that sounded like that, it's just the lowest note that we see in the Fourier transform is a peak corresponding to that low note and then all the overtones at about twice that frequency, there is another peak at about three times that frequency, there is an even smaller peak and all those overtones characterize the um.

The texture of a piece of music is what distinguishes, let's say, a human voice singing a note from a trumpet singing a note and all that and we can see it very well in the plot and in this case we don't have to take out only the first thousand terms we could do it across the entire signal, it's totally fine to plot that, but let's go ahead and limit our x values ​​again to say, let's look at those that go from zero to a thousand, I should say. zero to a thousand one to a thousand I'll do one to a thousand I just haven't put enough parentheses in, okay, now let me add to this graph the fast Fourier transform of the second file, the one that contains a note that is a little bit higher, that's why I'm going to make the first argument of this graph function be a list of different values, so I'm going to make it a list of the signal f of t yes, signal for t just say I because I in you know a range from one to two like this so let me see if I have my parentheses right here I think it should work and it's yelling at me no because it has an unexpected closing bracket oh yeah okay put it in too many what am I doing here?

I'm transmitting the absolute value through the ffts signal for i or for i in great, so here the blue graph gives us the Fourier transform of the first sound file that reads in that lowest note and then in the second one, you'll notice which has a peak above a slightly higher frequency and then all the harmonics for that which is nice and then if I do it for the first three sound files that we have then we see yet another frequency which is I think it was a musical sixth above the first and then all of its harmonics and then if we have the first four, uh, the first four files, then we also see the highest note appear.

Now the fun thing is if we actually took this plot and we were going to apply it to the sound file that contains all of them together, so instead of plotting four separate signals, we're just going to plot a single one, which is the last file that was. the fifth one on our list I still have to have my coma, you see all those same peaks but of course it doesn't know how to divide them in terms of colors because it doesn't know what files they come from, it just knows. that the one sound file you got, which sounds like they were all together, you let me go ahead and play them all together, so that's what you get and if you look at it just as a waveform, it's a pretty complicated thing, right?

You know? They are all these different waves added on top of each other, it's not entirely clear how to divide them into the different musical notes that make them up, but when you apply this Fourier transform to it you see these nice clean peaks on top of the relevant notes. and all the nuances of those notes, so with all that being said, let's go ahead and see how this discrete Fourier transform is mathematically defined, so I'm going to go ahead and show an animation here and the setup is that we're going to have some kind of signal which is a list of values, okay and for the moment we're just going to think aboutsomething really small, just eight different values ​​and here I'm plotting just a random sequence of values ​​and let's think about this as a list, I'll call the list s and at least in this context I think it's more natural to think of it as zero-indexed, so the initial element is at index zero, the next one is index one, index two, so and so on, so what I'm going to do is define for you what is the discrete Fourier transform of a signal like this , one term at a time and the first term is the easiest by far because it's just the sum of all these numbers and if we wanted to visualize it, you could think about taking all these vectors and just putting them together end to end and asking how long the result is . that's element zero of what I'm going to call shit and this is going to denote the Fourier transform, more specifically the discrete Fourier transform of the original signal, so it's pretty easy and things start to get a little more complex as we go.

I got to the other terms quite literally because we're going to start to visualize it on the complex plane and this seems like a lot at first, but we're going to break it down piece by piece. I have this Greek letter. zeta which I am using as a substitute for the term e to the negative pie 2 i divided by n now, whenever you see e raised to a purely imaginary number, this refers to a complex number that lies somewhere on the unit circle somewhere a distance one from the origin and what's in that exponent describes the angle, so in this context the capital n refers to the length of the list, so in this context n would be eight because there are eight elements in our list, so the angle would be negative two pi divided by eight, which is an eighth of a turn around the circle but walking clockwise, so is this vector sitting here an eighth of a turn and that would be zeta a the power one zeta to the zero power is the number one anything to the zero power is one and if we look at all the different powers of zeta, it gives us a good way to describe equally spaced angles around the unit circle and the idea of ​​transforms Fourier is to take some kind of signal and sort of walk around a circle but at several different frequencies, as you'll see at one point, and the language of complex numbers is introduced simply because that's the most mathematically natural way to describe Walking around a circle and doing all sorts of analysis makes getting through it so much easier, even if it initially seems intimidating to some to know that complex numbers have come into the game when it didn't seem like it should be part of reading in a sound file. , so all these coefficients, all these different zeta powers actually just say let's take our signal terms and let them rotate around a circle and our signal terms tell us how much we're going to scale these different vectors, like this which, for example, if we take, you know, The first element of our list multiplied by zeta raised to the zero power gives us a vector that is scaled appropriately and points just one unit to the right.

The next term gives us something that points in a southeast direction and so on walking around the circle exactly once. Now what we're doing is adding all of these different complex numbers, which could be thought of as adding a bunch of different vectors, and a good interpretation of that is to think about what the average value of all of these vectors is. um center of mass location for where all your tips are and then the sum as a whole will just be a scaled version of that and that can help you think about under what circumstances this total sum will be small versus when it's large. for example let me show you a circumstance where that sum would be quite large, let's say the original signal we took was something that looked like a cosine wave, it starts high and then halfway through the signal it becomes quite low and then as gets closer to the end of the sign, ends up high again, what that would mean for our vectors walking around the circle is that the one pointing to the right is quite long and then as you walk around the circle it gets more and more smaller. reaching a minimum by the time you get to the middle of the circle that would correspond to this fourth element here and then it starts to get bigger and bigger now when you take the sum of all of these because it's so asymmetrical that you're going to get something that's very large. and points to the right and the fact that it is very large would be telling you that this signal behaves like a cosine wave whose period is the entire length of the signal, starts high, goes down only once and ends high again. in that sense, this would be taking out a kind of low frequency term for the signal as a whole now if you took something else that has that same period, for example, something like a sine wave that starts at the mean value. goes up, then down, then up again, you get something that's twisted, it would be extracting that frequency, but if you think about what's happening with all the vectors, the point where they peak would occur at a different angle. then the resulting complex number that is obtained by adding all of these would again be very large, it would be far from the origin but it would have a different angle and what corresponds is the idea that this low frequency term that appears in the signal moves, has a different phase, so this complex number output it is important that it has two dimensions, a length and an angle, because that length tells us to what extent this signal contains the relevant frequency and the angle tells us what it is. the phase of that frequency as it appears, but those are just two terms of the discrete Fourier transform, let's start to look at the pattern by introducing other different signals and what they might produce there.

Okay, so with the second term we are doing something similar. but basically we're going to walk faster around the circle and this may seem a little confusing because the vectors are on top of each other, but let's go through one by one, the first term is the same as before. before it will always be this zeta to the power of 0, which is just too fancy a way of saying the number one, but we write it that way to fit the pattern, the next term here will be the square of that value, so basically, we're walking twice as fast as before, so instead of increasing our angle by negative 45 degrees, we increase it by negative 90 degrees and then the next term will be the fourth power, so here we are walking up. an even number zero two four six and so on and the fourth power in this context takes us all the way to the negative direction one, that's only because we have a small signal and there are only eight terms in it, so four is all that is you need to get halfway and generally you keep walking and this second term of the discrete Fourier transform corresponds to our vectors walking around the circle exactly twice, so again I want you to think about under what circumstances the sum of all these vectors would be that have now rotated in a slightly different way end up quite large and the quintessential example of a signal where that sum would be large would be again if we had a cosine wave, but this time one that has a higher frequency, so here would be if our Samples started high, went down and then went back up when you were in the middle of the signal and then went through the same cycle until going down and then going up to high because what that would mean is that we're walking around there.

Just so you know, up to a quarter of the way through the sign we already walk to the middle of the circle and then when we get to the middle of the circle it guides us completely, so if this was the behavior we would end up again. with a sort of unbalanced sum and again the fact that that sum is very large would be a reflection that this Fourier transform has captured a sort of pure frequency component within the original signal, so in this context, if we run this operation on a signal and we saw that the term at index 2 was really high, which would be telling us that there is a frequency within the signal that corresponds to two cycles over the duration of the signal and the pattern here essentially continues instead of having our powers from this basic complex.

The number increases one at a time or two at a time, increases by different numbers at a time and that will correspond to all the different elements within our Fourier transform, so the general formula seems to take a sum of all of them here. We're taking a sum where we're scaling the values ​​of the signal sn by some kind of power of this zeta term that I'm just using to more compactly write the e to the negative 2 pi i over n idea and all the thought. is that for different frequencies that you are trying to pick up, you have to walk faster around the circle, so for example, if it were a frequency of three, as you go through the sum, you would think that you walked three times around that circle and and so on, when the frequency is exactly half the length of the signal, we simply look back and forth between the points that are on the right side and on the left side and at five it walks even faster and again I know.

I'm being a little redundant, but each of these sums is a way of probing at a particular frequency that may or may not be part of the signal and as we saw by doing this with a really large signal that has a very clear frequency. sitting inside it because it was a musical note, you see a peak around that relevant frequency, so now that we have some kind of picture of what it's trying to do, let's see if we can implement this ourselves, let's write a transform function simple discrete Fourier and then we'll take a look at its runtime and see how it compares to this fast Fourier transform that comes to us from a library, so we'll define a function, let's just call it dft, that will receive a signal that we'll think of as a list of values, we'll think of them as real number values, but in principle they could also be complex number values, so we'll make sure to write it in a way that can take that into account.

Now it could begin. By defining this zeta value that we had, which exponentialized negative 2 pi i over n, that might be useful for us, so we'll say that zeta is the exponential of negative 2 pi. I guess I can write pi very well 2 pi multiplied by i, which in julia we can write as I am 2 pi i all divided by n and n in this context would be the length of the signal, okay, that gives us a zeta value now the end result we want will be some kind of list on the right and each element of that list looks like a sum, so the way I will do this is by describing the entire array at once using a list comprehension or array comprehension, each element It's going to look like some kind of sum and I'm going to go over the values ​​of f the frequency that go from 0 to n minus 1.

So the sum I'm also going to define using a type of understanding where it's going to be something that we're going to cover. the values ​​of n that goes from zero to n minus one and what we're adding that you can basically read from here is the nth element of our signal multiplied by a certain zeta power, so we can go ahead and say it's signal in n, but because in Julia things are indexed, I'm going to go ahead and add one to that and I guess I should use square brackets because it's not a function call.

I'm indexing into signal and then I'm going to multiply this by zeta. to the power of n times f Now, in principle, this should work, so let's go ahead and see if that actually ends up being the case. Invalid iteration specification, so for n in this range for f ah, I did the same thing again. Don't type the word, okay, now this feature is not optimal for several reasons that I'll point out in a moment, but we'll get to that in a moment, let's see if this actually does what we think. so let's say we're going to define a value signal that will be the first of our signals and let's see what the signal signal discrete Fourier transform does to it and it thinks and thinks and we can see that it's being very slow and it's not just slow because this is, You know, more than just the fast Fourier transform, there are very inefficient things that I've done with this initial attempt at implementation, so we can try to fix that in a moment, but in the meantime, let's go ahead and say, let's try with a signal smaller, something so let's just extract the first thousand terms from that and see how it works and unless I totally froze things with my bad implementation earlier, we should be able to. to get this working again, it looks like it's working and what we want to know is that this is doing the same thing as the function defined by the library.

Have we made some kind of mistake? If we just look at them, you know, it seems. like the first term is pretty similar, the second term is pretty similar and we could ask you if they're actually equal, so let's call this signal dft, we'll call it fft and we'll ask if they're equal to each other. andyou know, actually calculate them and it tells us no, they're not equal, but you know there's room for some numerical precision issues here and if we were to look at the difference between them, you know it looks like the first one is something that implies like a 10 to the power of 18 negative 10 to the power of 17 negative in general, you know the average value of the difference between these elements oh, we don't have the average, okay?

Let's say we're using statistics, let's use statistics, let's calculate the mean um and do the absolute value like this, it says I don't want to do the absolute value of a matrix because I have to convey it very well, so on average the difference is quite small, like this that they are basically the same and there is even a function we can call that gets more directly to what we want, namely approximately, where we can give it two different arrays and say, "Hey, are these two approximately the same?" and tells us that they are all. so our implementation has at least given us the right answer, but there are a couple of things that are not super efficient, so if we look at this runtime where you know we have it in an array that only has a length of 1000. our implementation , uh, Pluto tells us that it took 73 milliseconds and the speed The Fourier transform took 71 microseconds, so it's three orders of magnitude faster and if we increase the size, let's say it was something with a length of ten thousand instead of thousand, so it's ten times bigger, you can see it's actually taking quite a while, we have to wait. because here, you know, we're spending human time just on this simple discrete Fourier transform and it took eight seconds to implement, whereas the fast Fourier transform is still in microseconds, okay, so part of that is unnecessary inefficiencies that are not it's necessarily algorithmic where, for example, this zeta raised to a power is done quite uh n squared different times, you know, because we're doing a loop where we go from zero to n minus one and then another loop where we're doing the same thing, so if we have ten thousand correct terms, the square of that is a hundred million times.

We're doing something, so we really want to reduce what's happening 100 million times and in this case it's a little bit redundant. taking the zeta powers, the powers of our root of unity so many times because most of them are basically the same, there are only n different powers that could be because it involves walking around the circle with these uh n different steps, so what could do? Instead, what we do is we say, we're going to precompute a matrix of all the zeta powers where this is going to be e to the negative term 2 pi i, but we add a value n in there, we say that n goes from 0 to n. minus 1. and then what I want to be able to do is just say instead of calculating a power knowing that the CPU cycle runs there, I just want to look up the precomputed power, so I'll say zeta powers indexed at n times f but modulo n because there may be some loops that continue after we take it to the nth power and it goes back to one, except here I'm treating it as if it were a zero-indexed array when it's not because I want to be able to say things like zeta to the zeroth power should be the first element there, so what I'm going to do is actually convert this to an offset array, so that even though Julia is an indexed one, I can take a moment to pretend that it's zero-indexed at a time when that's a little bit more mathematically elegant and I basically just tell it what I want the indices of this to be, so I have the same array, but I tell it to make its indices go from 0 to n minus 1. and if that's how I define the things, then at least you don't do this exponential operation 100 million different times, and in fact, when we do that and run our thing again, previously we remember that it was eight seconds. run now is 322 milliseconds, much faster, but still nowhere near the fast Fourier transform that's still in those microseconds and you really start to feel this as soon as you scale things up because here you know we're going up to 10,000 elements. but if we go all the way and you know what the length of this array is, the length of the signal is, I guess we have to wait for it to compute the other things we're waiting for before we get the length is around 90,000, I think another increase of about an order of magnitude here we're seeing it again, it's human time and we have to wait, we're thinking in terms of seconds for this discrete Fourier transform to actually be calculated and this is even after I did oh, wait a minute, how is this working?

I just left this blank, zeta equals up there, did you really redefine the function the way you wanted? Now I'm wondering, well, in any case, I have to wait for everything. to load because this is the non-fast Fourier transform applied to a matrix of significant size, you know, 90,000 elements is not the largest you will find in the world, but it is of significant size and we are waiting for it, we are waiting for it, OK? So the fast Fourier transform took 14 milliseconds, so at least it got out of that microsecond range, but it's still a very fast operation and we're expecting just human seconds to pass twiddling our thumbs like the discrete transform.

Implemented the naive way takes its sweet time, so how long did that take? 42 seconds 42 seconds we had to sit here twiddling our thumbs while I computed the discrete transformation of a matrix with length 90,000 and this roughly makes sense again because we're doing these two different loops, each about something that's the length of the matrix, so if we take the square of that length, it means that this is the number of steps that need to be taken, so what is it that we are in millions when the 8 billion, so our implementation is either n squared and is the first one you might think of doing now.

What is a wonderful and somewhat surprising fact is that it is possible to do this discrete Fourier transform with an algorithm that is just o n log n so it scales much better. This is what is known as the fast Fourier transform. It's used so often that throughout computer science people don't even say discrete Fourier transform, they just say fft to refer to the mathematical operation because it's like why would you use any other algorithm? Ffts are Fourier transforms, and while I don't have time to go into detail about how that fft algorithm works here, I stumbled upon a video from a relatively young channel earlier this

week

. that gives a really beautiful overview of fast

fourier

transforms as an algorithm that gives you an idea of ​​where they come from, so it's a 30 minute video that's on a channel called reducible and I highly recommend it if you want to go into more detail. . details and one thing that's really nice is that he describes it in terms of polynomial multiplication and he motivates it with that and I don't know if you remember, but in one of the first lectures in this class I was talking about convolutions.

I mentioned what convolutions themselves are like, where you know you're imagining the signal traveling its way through an image and you're making this product of all the terms aligned with each other and adding them together that you can think of. which in terms of multiplying two polynomials, if you take one polynomial multiplied by another and ask what all the relevant terms are, a good way to think about it is to imagine reversing that first polynomial and somehow letting it cross over where the first one is. constant term just looks like multiplying the two constant terms, the next linear term involves multiplying two things, multiplying two things and adding them, and if you expand, you know the full product of these two polynomials and you really think about what's going on with all the different terms versus different powers of x, it is effectively a convolution between the two signals and it may not seem like a total coincidence that Fourier transforms are relevant here because if we look at the definition of our discrete Fourier transform, especially when we write with just zeta compacting all of this e to the two negative pi, each of these terms is essentially a polynomial um, it's a polynomial in zeta and we're evaluating it on several different possible inputs there, so to take a Fourier transform and evaluate a polynomial are essentially the same things as long as the points you're evaluating are certain specific roots of unity and it turns out you can take that as an idea to come up with this fast Fourier transform algorithm and again this video does it much better than...

I'll have time to do it here, so I'll go ahead and point that out to you, but it's a wonderful way to tie everything together.