Okay, so now I will show you how to implement, A discrete Fourier transform in SystemML. I have a file here, Which I'm saving to my home directory, called FT.dml, which will be a SystemML script. And I will start by first writing the gen_wave function that I used so far to show you different examples. Because I would like to have the same function in SystemML so I can generate test data. So I will give it the same name, gen_wave. Gen_wave is the function that accepts frequency and amplitude. It accepts a period, and also something great, which I will name Hz. Okay, I won't bother with phase shift, I just want to keep everything simple. And this function returns a matrix of doubles, called x. Okay, first I have to generate my x axis, my time. So I will generate the sequence that goes from 0 to t, and the resolution will be, the step will be t of the sampling created actually because I'm starting from 0, I will subtract 1 step. Okay, and then for x, I will take the amplitude, multiply it by sine, we remember the formula. Sine of 2pi. We don't actually have the pi constant in SystemML, so I will just pass it as an argument here. This will be a double. Okay, 2 times pi times the frequency times time. And again, I will just skip the phase shift. Now we have to define pi, which I will grab from my calculator. Where is pi? Pi is here. [SOUND] Okay, place it here. And then I will call the function, make sure that it does what I expect it to do. So let's say frequency of 2, amplitude of 1, a period of 150 Hertz and plus the pi value. Okay, and I want to test it now quickly. No, I don't need the calculator anymore. Okay, here's my SystemML folder. I will use the standalone version of SystemML. And now we'll pass my script as an argument. So this will be ft.dml. Okay, the script works, I don't see any output so maybe I'll just also print my matrix x. It's a matrix, so I have to convert it to string, this will be x, I will use two tops as separators. I don't need it in this case, but anyway. And I said precision, okay, and maybe I'll also get n which is a number of rows in x, and print this as well. Okay, so let me run the script again, what did I do? What did I do? It would help if I can type, decimal, okay. So I have 50 entries in x because my sampling rate is 50 over the same period. And then if you look at that the datapoints, you can actually see the wave. It starts 0, 0, I shift to 0, and then it goes all the way up to 1 and then it goes down to -1 and then it goes back to 0. Right, and I have two of those because my frequencies too. Now, I need to think about my direct Fourier transform implementation. So let me go back to the slides, right. And what I want to do now is implement these two sums, right, to give me. I'm not, Yeah. I need to implement these two sums here, this one and this one, to give me the Ak and Bk coefficients. And there are two ways of doing this. Probably the most simple one is if I look here, these k and n variables. I could actually create a loop over k, and then a loop over n, then compute Ak and Bk for each x of f, right. This is kind of, it's ugly, and it's not efficient. And we're using SystemML, so we would like to use metrics operations as much as possible because they'll be more efficient, and also more elegant. So, I will do something else. If I look at my algorithm here, I am missing an xn, [LAUGH] sorry. If I look at my algorithm here, I notice that these guys here, the argument of the sine function and the cosine function as well, don't actually depend on xn. If I get these indicator variables and create matrix, which I then multiply by two pi over n, I'm pretty much done. So what I can do is I can define a sequence, let's say n is a sequence that goes from 0 to N- 1, and K is also sequence that goes from 0 to N- 1, and these are vectors, so if I take the product of N times k transpose, this will give me a matrix of size N- 1 by N- 1. And this matrix, if I multiply by 2 pi over n, right, this will give me exactly this thing here. So then if I name this this matrix M, let's say. All I have to do is for my, Xa coefficient, And Xb coefficients, I can simply compute them as cosine of M, the product of the cosine of M and x. Right, and for Xb, the sine of M and x, which is my dataset, right? And this is all vectorized, so it will be more efficient and neater to implement. So let's go back to the script. Okay, and I will get rid of these two. Okay, and what I do now is I will define my two sequences from 0 to N- 1 step 1 and K is a sequence from 0 to N-1 step 1. And then my matrix M would be N times k transpose times 2 pi over n. Oops, right. And then Xa will be the cosine of M times x, and xp will be the product of sine M and, X. And I can merge these two together. I will create a matrix called DFT, and this will be Xa and Xb bonded together. And I can print this, just to see what the result is. So I'm printing the DFT matrix. Okay, let me, did it go, no, not you. I need my terminal which is here, okay? So let me run the script now. Okay, so this is my spectrogram in a way. I have a sampling rate of 50, so I get 50 bars. The frequency of the signal is 2, so that's 0, 1, 2, that's the second bar. And that's where I have my coefficients, and because the spectrogram is symmetrical, I also get this 25 there with the negative sign. Now, You might be wondering, why is this 25 here not equal to the amplitude, which is 1? And the reason is, because these are the Fourier coefficients, they don't reflect, they are not the same as the amplitude or the phase shift of the signal. So these are the coefficients that we use here, Find the correct slide, this one. Where did you go, here. These are the coefficients that we use here to actually reconstruct the signal, right? These are the Fourier coefficients. If we want to compute the magnitude of the signal and the phase shift, we have to do something else using these two coefficients again. So, if you want to compute the magnitude of the signal, let's say, For Xk, the magnitude for this frequency would be the Ak coefficient to the power of 2 plus the Bk coefficient to the power of 2, okay? And we take the square root of this, right, and this gives us the magnitude. For the phase shift I'll use capital fi. We need to take the inverse tangent of Bk over Ak, and this will give us the phase shift. So this is how we can compute this if we want to actually plot the signal or reconstruct the signal in a different way. So let's see if we can actually implement the magnitude in our script, and we have to remember that because this spectrogram, the Fourier coefficients are symmetric, We only care about the first half of the DFT matrix, right, because the other one is a mirror. So what we will have to do is we will get rid of the right hand side, we'll get rid of the trailing 25 elements in this case because my sample rate is 50. However, they still contain information about the magnitude, so we have to actually add them to the first 25 or just multiply the first 25 by 2 because there is symmetry, so it's the same thing. And then because we have sampling rate 50, we have to also divide by this sampling rate by number of data points to get the correct magnitude. So I will go back to my script here, and I will actually do this. So I will define a new matrix, let's say magnitude. This matrix will be the square root, but let me first take DFT, the first coefficient to the power of 2 plus the second coefficient to the power of 2, and the square root of this and then I will just take the half of this because I just need the first 25 elements. So from 1 to N over 2 and then I will multiply this by two and divide by N. Then I will print this matrix. Okay, where is my terminal? Let's run through the script again. I made a mistake. DFT, okay. Okay, so now I have the magnitude. And as you can see, it's absolutely spot on. Does the third bin which corresponds to frequency of 2 hertz and a magnitude of 1. And that's just the left-hand side of the spectrogram, okay? So this was our, Discrete Fourier transformation implementation in SystemML.
