Fourier Transform
Introduction
The Fourier Transform is an important tool in Image Processing, and is directly related to filter theory, since a filter, which is a convolution in the spatial domain (=the image), is a simple multiplication in the spectral domain (= the FT of the image)!
Most other tutorials about Fourier Transforms of images are in boring greyscale. Even, many programs and plugins available that can calculate the FT on images only support greyscale. In times where 24 and 32 bit color are getting old, we cannot accept this any longer! This tutorial is in 24-bit color, by doing the FT on each colorchannel separately. Here's an example where the FT of an image is used to make a tillable grass texture look better:
There is also another way to calculate the FT of color images however, by using quaternions instead of complex numbers, which can convert an image of 4 channels to another image with 4 channels.
Before beginning with the Fourier Transform on images, which is the 2D version of the FT, we'll start with the easier 1D FT, which is often used for audio and electromagnetical signals.
To fully understand this tutorial, it'd be interesting to know about complex numbers, their magnitude (amplitude), and argument (phase). If you'd like to learn more about these, there's an appendix about complex numbers (one of the other html files), or find a tutorial with google.
Signals and Spectra
Before starting with the Fourier Transform, an introduction about Signals and Spectra is in place. As signals, for now consider signals in the time, like audio and electromagnetic signals. If you plot such a signal in function of time, a simple signal could look like this:
This is a sine wave, so if this'd be an audio signal, you'd hear a pure tone with only 1 frequency (like the beep of a PC speaker, or a very pure flute). The grey curve behind the red curve represents the amplitude of the signal, which is always positive.
The spectrum of a signal contains for every frequency, how much of this frequency is in the signal. Since the signal above is a sine, it has only one frequency, so the spectrum would be very simple: only one value will be positive (= the frequency of the sine curve), and all others are zero. So the spectrum would have a single peak. A spectrum has two sides however: the negative side on the left, and the positive side on the right. The negative side contains negative frequencies. For real signals (with no imaginary part), like audio signals are, the negative side of the spectrum is always a mirrored version of the positive side. So for the sine signal above, the positive side will have a single peak, and this peak will be mirrored in the negative side. This is what it's spectrum looks like:
The white curve represents the amplitude of the signal. The red and green curve represent the real and imaginary part of the spectrum, respectively, but a spectrum is raraly represented that way. Normally, both the amplitude and phase of a spectrum are given, but since there's not that much interesting info for us in it now, the phase isn't included.
Let's study another example now: the same sine, but with a DC component added: sin(x)+C. The name DC component comes from electronics, where it stands for Direct Current (as opposed to Alternating Current). The DC component in any signal, is the time average of it. A sine signal has time average 0, but if however, you add a constant to the sine, the signals time average will become that constant. This adding of the constant, results in the sine signal being raised higher above the horizontal axis:
Such a constant, or DC component, has frequency 0. So we can expect that the spectrum will get an extra peak at the origin, and indeed, it has:
So, if the spectrum of a signal is non-zero at the origin, you know that it's time average isn't zero and it has a DC component! The further away from the origin to the left or right a peak is, the higher the frequency it represents, so a peak far to the right (and left) means that the signal contains a high frequency component. Note that an audio signal never has a DC component, since nobody can produce or hear it. Electrical signals can have it however.
Now let's look at a signal which is the sum of two sine functions: the second sine function has the double frequency of the first, so the curve has the shape sin(x)+sin(2*x):
Since there are now two sines with two different frequencies, we can expect two peaks on the positive side of the spectrum (and two more on the negative side since it's the mirrored version):
If the two sines both have a different phase (i.e. one of the two sines is shifted horizontally), the amplitude of the spectum will still be the same, only the phase will be different, but the phase isn't shown here.
Arbitrary sound signals, such as a spoken word, consist of an infinite sum (or integral) of infinitely much sine functions, all with a different frequency, and all with their own amplitude and phase, and the spectrum is the perfect plot to view how much of each frequency is in the signal. Such a spectrum won't just have a few peaks, but will be a continuous function.
There are a few special functions with a special spectrum:
The Dirac impuls is a signal that is zero everywhere, except at the origin, where it's infinite. That is the ideal version for continuous functions, for a discrete function on a computer it can as well be represented as a single peak with finite height at the origin.
Such a peak sounds like a popping sound in an audio signal, and contains ALL frequencies. That's why it's spectrum looks like this (the horizontal black line):
The spectrum is positive everywhere, so every frequency is contained in the signal. This also means that, to get such a peak as a sum of sine functions, you have to add infinitely much elementary sine functions all with the same amplitude, and all shifted with a certain phase, they'll cancel each other out, except at the origin, hence the peak.
Note that it's possible to intepret the spectrum as a time signal (a constant DC signal), and that then the red function can be seen as it's spectrum (a peak at zero, since a DC signal has frequency 0). This duality is one of the interesting properties of the Fourier Transform as we'll see later.
Finally, another special signal is the sinc(x) function; sinc(x) = sin(x)/x:
It's spectrum is a rectangular pulse (here only an approximation of it because the spectrum was calculated using only 128 points):
Again because of the duality between a signal and it's spectrum, a rectangular time signal will have a sinc function as it's spectrum.
Spectra are also used to study the color of light. These are exactly the same thing, since light too is a 1-dimensional signal containing different frequencies, some of which the eye happens to be sensitive to. The spectrum shows for every frequency, how much of it is in the light, just the same as the audio spectrum does for audio signals.
An MP3 player like Winamp also shows the spectrum of the audio signal it's playing (but it varies with time because they calculate the spectrum each time again for the last few milliseconds), and spectra of audio signals can be used to detect bass, treble and vocals in music, for speach recognition, audio filters, etc... Since this is a tutorial about computer graphics, I won't go any deeper in this however.
Often, the spectrum is drawn on a logarithmical scale instead of a linear one like here, then the amplitude is expressed in decibels (dB).
The Fourier Transform
So how do you calculate the spectrum of a given signal? With the Fourier Transform!
The Continuous Fourier Transform, for use on continuous signals, is defined as follows:
And the Inverse Continuous Fourier Transform, which allows you to go from the spectrum back to the signal, is defined as:
F(w) is the spectrum, where w represents the frequency, and f(x) is the signal in the time where x represents the time. i is sqrt(-1), see complex number theory.
Note the similarity between both transforms, which explains the duality between a signal and it's spectrum.
A computer can't work with continuous signals, only with finite discrete signals. Those are finite in time, and have only a discrete set of points (i.e. the signal is sampled). One of the properties of the Fourier Transform and it's inverse is, that the FT of a discrete signal is periodic. Since for a computer both the signal and the specrum must be discrete, both the signal and spectrum will be periodic. But by taking only one period of it, we get our finite signal. So if you're taking the DFT of a signal or an image on the computer, mathematicaly speaking, that signal is infinitely repeated, or the image infinitely tiled, and so is the spectrum. A nice property is that both the signal and the spectrum will have the same number of discrete points, so the DFT of an image of 128*128 pixels will also have 128*128 pixels.
Since the signal is finite in time, the infinite borders of the integrals can be replaced by finite ones, and the integral symbol can be replaced by a sum. So the DFT is defined as:
And the inverse as:
This looks already much more programmable on a computer, to program the FT, you have to do the calculation for every n, so you get a double for loop (one for every n, and one for every k, which is the sum). The exponential of the imaginary number can be replaced by a cos and a sin by using the famous formuma e^ix = cosx + isinx. More on programming it follows further in this tutorial.
There are multiple definitions of the DFT, for example you can also do the division through N in the forward DFT instead of the inverse one, or divide through sqrt(N) in both. To plot it on a computer, the best result is gotten by dividing through N in the forward DFT.
Properties of The Fourier Transform
A signal is often denoted with a small letter, and it's fourier transform or spectrum with a capital letter. The relatoin between a signal and it's spectrum is often denoted with f(x) <--> F(w), with the signal on the left and it's spectrum on the right.
The Fourier Transform has some interesting properties, some of which make it easier to understand why spectra of certain signals have a certain shape.
Don't take this as a proper list however, some scaling factors like 2*pi are left away, they're not that important and can be left away or depending if you use frequency or pulsation. Handbooks of Sigal Processing and sites like Wikipedia and Mathworld contain much more complete and mathematically correct tables of FT properties. We're most interested in the shape of the functions and not the scale.
Linearity
f(x)+g(x) <--> F(w)+G(w)
a*f(x) <--> a*F(w)
This means that if you add/substract two signals, their spectra are added/substracted as wall, and if you increase/decrease the amplitude of the signal, the amplitude of it's spectrum will be increased/decreased with the same factor.
Scaling
f(a*x) <--> (1/a) * F(w/a)
This means that if you make the function wider in the x-direction, it's spectrum will become smaller in the x-direction, and vica versa. The amplitude will also be changed.
Time Shifting
f(x-x0) <--> (exp(-i*w*x0)) * F(w)
Since the only thing that happens if you shift the time, is a multiplication of the Fourier Transform with the exponential of an imaginary number, you won't see the difference of a time shift in the amplitude of the spectrum, only in the phase.
Frequency Shifting
(exp(-i*w0*x))*f(x)<-->F(w-w0)
This is the dual of the time shifting.
Duality or Symmetry
if f(x) <--> F(w)
then F(x) <--> f(-w)
Apart from some scaling factors at least. Because of this property, for example, the spectrum of a rectangular pulse is a sinc function and at the same time the spectrum of a sinc function is a rectangular pulse.
Symmetry Rules
These are only a few of the symmetry rules:
- The Fourier Transform of a real even signal is real and even (i.e. it's symmetrical, mirrored around the y-axis)
- The Fourier Transform of a real odd signal is imaginary and odd (odd means it's assymetrical, mirrored around the origin point)
- Since arbitrary real signals are always a sum of an even and an odd function, the Fourier Transform of a real signal has an even real part and an odd imaginary part, and the Amplitude is thus always symmetrical.
- In the same way, the amplitude of FT's of pure imaginary signals are symmetrical, but the FT of complex signals isn't always symmetrical.
Convolution is an operation between two functions that is defined as an integral, but can also be explained as follows:
Take 2 functions, for example two rectangle functions f1(x) and f2(x). Keep one of the rectangles at a fixed position. Mirror the second one around the y axis (this doesn't have a lot of effect on a rectangle). Then shift this second function over a value u, where u goes from -infinity to +infinity.
The resulting function g(u), the convolution of f1(x) and f2(x), takes u as a parameter. For a certain u, multiply f1 and the shifted and mirrored f2 with each other, and take area under the result (by integrating it). This area is the value of g(u) for that u! This has to be done for every u to know the complete result. For example, for the two rectangle functions, as we start at u=-infinity, g(u) is 0 because the two rectangles don't overlap, and multiplication of the two functions will thus result in the zero function, which has zero area. As u moves to the right, g(u) will remain 0 all the time, until finally the two rectangles start overlapping. Now the more we go to the right, the larger the area will become, until it reaches it's maximum, so g(u) is rising. Then, the rectangles are overlapping less and less again as the second one moves more to the right, and the value of g(u) decreases again. Finally g(u) reaches zero again and remains so until +infinity. The resulting g(u) is thus a triangle function.
This is what filters of painting programs that don't use the FT do in 2D. The image is then the function that remains at fixed position, while the filter matrix moves around the whole image, and you calculate the value of every pixel by multiplying every element of the matrix with every corresponding pixel of the image the matrix overlaps. For large filter matrices this is a lot of work. You can create your own filter matrices in painting programs like Paint Shop Pro and Photoshop by using the "User Defined" filter.
Correlation is almost the same as convolution, except you don't have to mirror the second function.
The convolution theorem sais that, convolving two functions in the time domain corresponds to multiplying their spectra in the fourier domain, and vica versa. This multiplication is a much easier operation than convolution.
Since filters, too, can be described in the time and the fourier domain (called the impulse response and transfer function of the filter respectively), and linear filters do a convolutoin in the time domain, in the frequency domain this corresponds to multiplying the transfer function of the filter with the spectrum of the signal. The impulse response of the filter in the time domain is the response of the filter to a dirac impuls (see earlier), and the transfer function of the filter is the Fourier Transform of the impulse response.
For us, the Convolution Theorem will come in handy when we experiment with the Fourier Transformations of signals and images.
Programming The 1D Fourier Transform
Here's the code of a simple program that'll calculate the FT of a given signal, and will plot both the signal and the FT. Put it in the main.cpp file.
First the variables and functions are declared. N will be the number of discrete points in the signal. For the signal and it's spectrum, respectively a small f or g and capital F or G are used, as it's also done in mathematics.
const int N = 128; //yeah ok, we use old fashioned fixed size arrays here
double fRe[N]; //the function's real part, imaginary part, and amplitude
double fIm[N];
double fAmp[N];
double FRe[N]; //the FT's real part, imaginary part and amplitude
double FIm[N];
double FAmp[N];
const double pi = 3.1415926535897932384626433832795;
void DFT(int n, bool inverse, double *gRe, double *gIm, double *GRe, double *GIm);
//Calculates the DFT
void plot(int yPos, double *g, double scale, bool trans, ColorRGB color);
//plots a function
void calculateAmp(int n, double *ga, double *gRe, double *gIm);
//calculates the amplitude of a complex function
|
The main function creates a sine signal, and will then use the other functions to calculate the FT and plot it.
int main(int /*argc*/, char */*argv*/[])
{
screen(640, 480, 0, "1D DFT");
cls(RGB_White);
for(int x = 0; x < N; x++)
{
fRe[x] = 25 * sin(x / 2.0); //generate a sine signal. Put anything else you like here.
fIm[x] = 0;
}
//Plot the signal
calculateAmp(N, fAmp, fRe, fIm);
plot(100, fAmp, 1.0, 0, ColorRGB(160, 160, 160));
plot(100, fIm, 1.0, 0, ColorRGB(128, 255, 128));
plot(100, fRe, 1.0, 0, ColorRGB(255, 0, 0));
DFT(N, 0, fRe, fIm, FRe, FIm);
//Plot the FT of the signal
calculateAmp(N, FAmp, FRe, FIm);
plot(350, FRe, 12.0, 1, ColorRGB(255, 128, 128));
plot(350, FIm, 12.0, 1, ColorRGB(128, 255, 128));
plot(350, FAmp, 12.0, 1, ColorRGB(0, 0, 0));
drawLine(w / 2, 0, w / 2, h - 1, ColorRGB(128, 128, 255));
redraw();
sleep();
return 0;
}
|
Then comes the function that does the actual work (albeit with a very suboptimal algorithm here): the DFT function
void DFT(int n, bool inverse, double *gRe, double *gIm, double *GRe, double *GIm)
{
for(int w = 0; w < n; w++)
{
GRe[w] = GIm[w] = 0;
for(int x = 0; x < n; x++)
{
double a = -2 * pi * w * x / float(n);
if(inverse) a = -a;
double ca = cos(a);
double sa = sin(a);
GRe[w] += gRe[x] * ca - gIm[x] * sa;
GIm[w] += gRe[x] * sa + gIm[x] * ca;
}
if(!inverse)
{
GRe[w] /= n;
GIm[w] /= n;
}
}
}
|
It takes 4 arrays as parameters: *gRe and *gIm to read from (these are the real and imaginary part of the signal), and *GRe and *GIm to write the results to (these are the FT of the signal). Both the input and output are complex numbers, which is why two arrays (a real and an imaginary one) are needed. n is the size of the signal (and the arrays), so it's the number of discrete points in the signal.
The first for loop with w going from 0 to n-1 means that a calculation is needed for every single point of the result. The second for loop, with x going from 0 to n, means that every point of the input signal needs to be taken in account, it's the sum from the mathematical definition of the DFT. The exponential function is changed by the formula e^ix = cosx + isinx, and the real and imaginary part are calculated seperatly. Afterwards the values are divided through n, so that the values are reasonably small enough to plot them on the screen.
By setting "inverse" to true, it can also calculate the IDFT (Inverse Discrete Fourier Transform), but this isn't used in this example. The only differences are a different sign in the exponential function and no division through n.
Next, the plot function takes a single array as parameter, and will draw the signal or spectrum contained in this array on screen with a color you want using lines and filled circles. You can also use a parameter "shift" to bring the edges to the center and vica versa. This is useful for a better representation of the Fourier Transform: normally the result of the DFT has the low frequency components on the sides, and the highest frequency in the center. Since the result is in fact periodic, bringing the sides to the center is the same as moving the plot half a period.
void plot(int yPos, double *g, double scale, bool shift, ColorRGB color)
{
drawLine(0, yPos, w - 1, yPos, ColorRGB(128, 128, 255));
for(int x = 1; x < N; x++)
{
int x1, x2, y1, y2;
int x3, x4, y3, y4;
x1 = (x - 1) * 5;
//Get first endpoint, use the one half a period away if shift is used
if(!shift) y1 = -int(scale * g[x - 1]) + yPos;
else y1 = -int(scale * g[(x - 1 + N / 2) % N]) + yPos;
x2 = x * 5;
//Get next endpoint, use the one half a period away if shift is used
if(!shift) y2 = -int(scale * g[x]) + yPos;
else y2 = -int(scale * g[(x + N / 2) % N]) + yPos;
//Clip the line to the screen so we won't be drawing outside the screen
clipLine(x1, y1, x2, y2, x3, y3, x4, y4);
drawLine(x3, y3, x4, y4, color);
//Draw circles to show that our function is made out of discrete points
drawDisk(x4, y4, 2, color);
}
}
|
And finally, the function calculateAmp calculates the amplitude given the real and imaginary part, so you can plot the amplitude too.
//Calculates the amplitude of *gRe and *gIm and puts the result in *gAmp
void calculateAmp(int n, double *gAmp, double *gRe, double *gIm)
{
for(int x = 0; x < n; x++)
{
gAmp[x] = sqrt(gRe[x] * gRe[x] + gIm[x] * gIm[x]);
}
}
|
The output of this program is as follows:
![](https://i-blog.csdnimg.cn/blog_migrate/8a004132915ada150e87177431953343.gif)
The red function is the input function, which is a sine, and the black funtion at the bottom is the amplitude of it's spectrum, which is calculated by the DFT function. The green curves are the imaginary parts.
The Fast Fourier Transform
The above DFT function correctly calculates the Discrete Fourier Transform, but uses two for loops of n times, so it takes O(n�) arithmetical operations. A faster algorithm is the Fast Fourier Transform or FFT, which uses only O(n*logn) operations. This makes a big difference for very large n: if n would be 1024, the DFT function would take 1048576 (about 1 million) loops, while the FFT would use only 10240.
The FFT works by splitting the signal into two halves: one halve contains all the values with even index, the other all the ones with odd index. Then, it splits up these two halves again, etc..., recursivily. Due to this, a requirement is that the size of the signal, n, must be a power of 2. There exist other FFT versions that can work on signals that can have length n*m where both n and m are primes, but that goes beyond the scope of this tutorial.
Here's a new version of the DFT function, called FFT, which uses the Fast Fourier Transform instead. The FFT algorithm used here is called the Radix-2 algorithm, by Cooley-Tukey, developed 1965. You can add it to the previous program and change the function call of it to FFT to test it out.
First, this function will calculate the logarithm of n, needed for the algorithm.
void FFT(int n, bool inverse, double *gRe, double *gIm, double *GRe, double *GIm)
{
//Calculate m=log_2(n)
int m = 0;
int p = 1;
while(p < n)
{
p *= 2;
m++;
}
|
Next, Bit Reversal is performed. Because of the way the FFT works, by splitting the signal in it's even and odd components every time, the algorithm requires that the input signal is already in the form where all odd components come first, then the even ones, and so on in each of the halves. This is the same as reversing the bits of the index of each value (for example the second value, with index 0001 will now get index 1000 and this go to the center), hence the name.
//Bit reversal
GRe[n - 1] = gRe[n - 1];
GIm[n - 1] = gIm[n - 1];
int j = 0;
for(int i = 0; i < n - 1; i++)
{
GRe[i] = gRe[j];
GIm[i] = gIm[j];
int k = n / 2;
while(k <= j)
{
j -= k;
k /= 2;
}
j += k;
}
|
Next the actual FFT is performed. Mathematically speaking it's a recursive function, but this is a modified version to work non-recursively, it'd be way too hard to pass complete arrays to a recursive function every time. ca and sa still represent the cosine and sine, but are calculated with the half-angle formulas and initual values -1.0 and 0.0 (cos and sin of pi): each loop, ca and sa represent the cosine and sine of half the angle as the previous loop.
//Calculate the FFT
double ca = -1.0;
double sa = 0.0;
int l1 = 1, l2 = 1;
for(int l = 0; l < m; l++)
{
l1 = l2;
l2 *= 2;
double u1 = 1.0;
double u2 = 0.0;
for(int j = 0; j < l1; j++)
{
for(int i = j; i < n; i += l2)
{
int i1 = i + l1;
double t1 = u1 * GRe[i1] - u2 * GIm[i1];
double t2 = u1 * GIm[i1] + u2 * GRe[i1];
GRe[i1] = GRe[i] - t1;
GIm[i1] = GIm[i] - t2;
GRe[i] += t1;
GIm[i] += t2;
}
double z = u1 * ca - u2 * sa;
u2 = u1 * sa + u2 * ca;
u1 = z;
}
sa = sqrt((1.0 - ca) / 2.0);
if(!inverse) sa =- sa;
ca = sqrt((1.0 + ca) / 2.0);
}
|
Finally, the values are divided through n if it isn't the inverse DFT.
//Divide through n if it isn't the IDFT
if(!inverse)
for(int i = 0; i < n; i++)
{
GRe[i] /= n;
GIm[i] /= n;
}
}
|
The function reads from *gRe and *gIm again and stores the result in *GRe and *GIm.
There are 3 nested loops now in the FFT part, but only the first one goes from 0 to n-1, the other two together only loop log_2(n) times. For this small signal of only 128 values you probably won't notice any increase in speed because both DFT and FFT are calculated instantly (unless you're working on some computer from the '40s), but on bigger signals, on 2D signals, and in applications where the FT has be be calculated over and over in real time, it'll make an important difference.
Filters
In the Fourier Domain, it's easy to apply some filters. All you have to do is multiply the spectrum with the transfer function of the filter.
Low pass filters only let through components with low frequencies (for example the bass from music), while High Pass filters let only components with high frequencies through. You can also make Band Pass and Band Stop filters, which let through or block components with a certain frequency. Amplifiers multiply the whole spectrum with a constant value, and you can of course make filters that amplify certain frequencies, or triangular filters that weaken higher frequencies more and more, etc...
So to make for example a Low Pass filter, multiply the spectrum with a rectangular function with the rectangle around the origin. All the low frequency components will then be kept, while the high ones are multiplied by zero and are thus filtered away.
Here's a part of code for a program that'll first allow you to select different signals and watch their FT, and after pressing the space key, will allow you to choose different filters on the spectrum, and it'll calculate the inverse FT to see how the original signal is affected by the filter.
First all variables and functions are declared.
const int N = 128; //yeah ok, we use old fashioned fixed size arrays here
double fRe[N]; //the function's real part, imaginary part, and amplitude
double fIm[N];
double fAmp[N];
double FRe[N]; //the FT's real part, imaginary part and amplitude
double FIm[N];
double FAmp[N];
const double pi = 3.1415926535897932384626433832795;
void FFT(int n, bool inverse, double *gRe, double *gIm, double *GRe, double *GIm); //Calculates the DFT
void plot(int yPos, double *g, double scale, bool trans, ColorRGB color); //plots a function
void calculateAmp(int n, double *ga, double *gRe, double *gIm); //calculates the amplitude of a
complex function
|
The main function starts here. It enters the first loop, this loop will allow you to choose a function. The parameters p1 and p2 can be changed to modify the amplitude and/or width of some of the functions.
int main(int /*argc*/, char */*argv*/[])
{
screen(640, 480, 0, "Fourier Transform and Filters");
for(int x = 0; x < N; x++) fRe[x] = fIm[x] = 0;
bool endloop = 0, changed = 1;
int n = N;
double p1 = 25.0, p2 = 2.0;
while(!endloop)
{
if(done()) end();
readKeys();
|
This is inside the loop after it has read the keys. Every key has another action assiciated to it to create a certain function, for example if you press the g key (SDLK_g) it'll put a sum of two sines in fRe[x]. The arrow keys also have some actions to set the imaginary part to 0, equal to the real part, or a modified version of the real part. Finally, the space keys sets "endloop" to true so that the first loop will end.
Not all keys are included here, because it's very messy and big code. Only a few are shown as an example, the rest is in the downloadable source file.
if(keyPressed(SDLK_a)) {for(int x = 0; x < n; x++) fRe[x] = 0; changed = 1;} //zero
if(keyPressed(SDLK_b)) {for(int x = 0; x < n; x++) fRe[x] = p1; changed = 1;} //constant
if(keyPressed(SDLK_c)) {for(int x = 0; x < n; x++) fRe[x] = p1 * sin(x / p2); changed = 1;} //sine
//ETCETERA... *snip*
if(keyPressed(SDLK_SPACE)) endloop = 1;
|
This is the last part of the first loop, where it draws the chosen function, calculates it's FT, and draws that one too.
if(changed)
{
cls(RGB_White);
calculateAmp(N, fAmp, fRe, fIm);
plot(100, fAmp, 1.0, 0, ColorRGB(160, 160, 160));
plot(100, fIm, 1.0, 0, ColorRGB(128, 255, 128));
plot(100, fRe, 1.0, 0, ColorRGB(255, 0, 0));
FFT(N, 0, fRe, fIm, FRe, FIm);
calculateAmp(N, FAmp, FRe, FIm);
plot(350, FRe, 12.0, 1, ColorRGB(255, 128, 128));
plot(350, FIm, 12.0, 1, ColorRGB(128, 255, 128));
plot(350, FAmp, 12.0, 1, ColorRGB(0, 0, 0));
drawLine(w / 2, 0, w / 2, h - 1, ColorRGB(128, 128, 255));
print("Press a-z to choose a function, press the arrows to change the imaginary part,
press space to accept.", 0, 0, RGB_Black);
redraw();
}
changed = 0;
}
|
After the first loop has ended, the second part of the main function starts. First, an extra spectrum buffer is created, so that both the original spectrum and it's filtered version can be remembered. Then the second loop starts.
//The original FRe2, FIm2 and FAmp2 will become multiplied by the filter
double FRe2[N];
double FIm2[N];
double FAmp2[N];
for(int x = 0; x < N; x++)
{
FRe2[x] = FRe[x];
FIm2[x] = FIm[x];
}
cls();
changed = 1; endloop = 0;
while(!endloop)
{
if(done()) end();
readKeys();
|
Now the keys all change the spectrum in a certain way: they'll aplly a filter. Low Pass filters (LP) set the spectrum to 0, except around the origin. High Pass filters (HP) do the oposite.
Not all keys are shown here, because it's very messy and big code. Only a few are shown as an example, the rest is in the downloadable source file in a more closely packed form.
if(keyPressed(SDLK_a))
{
for(int x = 0; x < n; x++)
{
FRe2[x] = FRe[x];
FIm2[x] = FIm[x];
}
changed = 1;
} //no filter
if(keyPressed(SDLK_b))
{
for(int x = 0; x < n; x++)
{
if(x < 44 || x > N - 44)
{
FRe2[x] = FRe[x];
FIm2[x] = FIm[x];
}
else FRe2[x] = FIm2[x] = 0;
}
changed = 1;
} //LP
//ETCETERA... *snip*
|
Again, if a key was pressed, the new spectrum is drawn, and the inverse FFT is performed to draw the filtered version of the original function. If you press escape the program will end.
if(changed)
{
cls(RGB_White);
FFT(N, 1, FRe2, FIm2, fRe, fIm);
calculateAmp(N, fAmp, fRe, fIm);
plot(100, fAmp, 1.0, 0, ColorRGB(160, 160, 160));
plot(100, fIm, 1.0, 0, ColorRGB(128, 255, 128));
plot(100, fRe, 1.0, 0, ColorRGB(255, 0, 0));
calculateAmp(N, FAmp2, FRe2, FIm2);
plot(350, FRe2, 12.0, 1, ColorRGB(255, 128, 128));
plot(350, FIm2, 12.0, 1, ColorRGB(128, 255, 128));
plot(350, FAmp2, 12.0, 1, ColorRGB(0, 0, 0));
drawLine(w/2, 0,w/2,h-1, ColorRGB(128, 128, 255));
print("Press a-z to choose a filter. Press esc to quit.", 0, 0, RGB_Black);
redraw();
}
changed=0;
}
return 0;
}
|
The functions FFT, plot and calculateAmp are identical to the ones explained earlier and aren't included here anymore.
There's not much new in this code, the functions FFT, plot and calculateAmp are still the same, but by playing with this program you can see the effects of Low Pass (LP), High Pass (HP) and a few other filters. Just follow the instructions on screen.
For example, here, a signal which is the sum of 5 sines and a DC component (gotten by pressing the "p" key) is shown together with it's spectrum, and then, respectively, a LP and a HP are applied to it:
The original signal:
![](https://i-blog.csdnimg.cn/blog_migrate/5be1a239aea9cd49ecd094e842b6b563.gif)
An LP Filter:
![](https://i-blog.csdnimg.cn/blog_migrate/bdb98151b7a36571f87833df225778f9.gif)
A HP Filter
![](https://i-blog.csdnimg.cn/blog_migrate/872fbd61c576982b7d0cfa696e01be32.gif)
Here, the same is done with a step function (gotten by pressing the "e" key):
![](https://i-blog.csdnimg.cn/blog_migrate/ce8f5e912f70878ab3fe745493388aa8.gif)
An LP filter on it:
![](https://i-blog.csdnimg.cn/blog_migrate/2edeb05812a2ae7c0be5171d8095b399.gif)
A HP filter on it. This filter filters out almost only the DC component, so two almost-flat lines remain in the signal:
![](https://i-blog.csdnimg.cn/blog_migrate/3e82094adf19ea9c21991fc3a835a34e.gif)