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

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;
}

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;
}
}
}

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);
}
}

//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]);
}
}

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++;
}

//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;
}

//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);
}

//Divide through n if it isn't the IDFT
if(!inverse)
for(int i = 0; i < n; i++)
{
GRe[i] /= n;
GIm[i] /= n;
}
}

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

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();

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;

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 az to choose a function, press the arrows to change the imaginary part, press space to accept.", 0, 0, RGB_Black);
redraw();
}
changed = 0;
}

//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();

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*

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,h1, ColorRGB(128, 128, 255));
print("Press az to choose a filter. Press esc to quit.", 0, 0, RGB_Black);
redraw();
}
changed=0;
}
return 0;
}

//yeah, we're working with fixed sizes again...
const int N = 128; //the width of the image
const int M = 128; //the height of the image
double fRe[N][M][3], fIm[N][M][3], fAmp[N][M][3]; //the signal's real part, imaginary part, and amplitude
double FRe[N][M][3], FIm[N][M][3], FAmp[N][M][3]; //the FT's real part, imaginary part and amplitude
double fRe2[N][M][3], fIm2[N][M][3], fAmp2[N][M][3]; //will become the signal again after IDFT of the spectrum
double FRe2[N][M][3], FIm2[N][M][3], FAmp2[N][M][3]; //filtered spectrum
double pi = 3.1415926535897932384626433832795;
void draw(int xpos, int yPos, int n, int m, double *g, bool shift, bool neg128);
void DFT2D(int n, int m, bool inverse, double *gRe, double *gIm, double *GRe, double *GIm);
void calculateAmp(int n, int m, double *ga, double *gRe, double *gIm);

int main(int /*argc*/, char */*argv*/[])
{
screen(N * 3, M * 3, 0, "2D DFT and FFT");
std::vector<ColorRGB> img;
unsigned long dummyw, dummyh;
if(loadImage(img, dummyw, dummyh, "pics/test.png"))
{
print("image pics/test.png not found");
redraw();
sleep();
cls();
}
//set signal to the image
for(int x = 0; x < N; x++)
for(int y = 0; y < M; y++)
{
fRe[x][y][0] = img[N * y + x].r;
fRe[x][y][1] = img[N * y + x].g;
fRe[x][y][2] = img[N * y + x].b;
}
//draw the image (real, imaginary and amplitude)
calculateAmp(N, M, fAmp[0][0], fRe[0][0], fIm[0][0]);
draw(0, 0, N, M, fRe[0][0], 0, 0);
draw(N, 0, N, M, fIm[0][0], 0, 0);
draw(2 * N, 0, N, M, fAmp[0][0], 0, 0);
//calculate and draw the FT of the image
DFT2D(N, M, 0, fRe[0][0], fIm[0][0], FRe[0][0], FIm[0][0]); //Feel free to change this to FFT2D
calculateAmp(N, M, FAmp[0][0], FRe[0][0], FIm[0][0]);
draw(0, M, N, M, FRe[0][0], 1, 1);
draw(N, M, N, M, FIm[0][0], 1, 1);
draw(2 * N, M, N, M, FAmp[0][0], 1, 0);
//calculate and draw the image again
DFT2D(N, M, 1, FRe[0][0], FIm[0][0], fRe2[0][0], fIm2[0][0]); //Feel free to change this to FFT2D
calculateAmp(N, M, fAmp2[0][0], fRe2[0][0], fIm2[0][0]);
draw(0, 2 * M, N, M, fRe2[0][0], 0, 0);
draw(N, 2 * M, N, M, fIm2[0][0], 0, 0);
draw(2 * N, 2 * M, N, M, fAmp2[0][0], 0, 0);
redraw();
sleep();
return 0;
}

void DFT2D(int n, int m, bool inverse, double *gRe, double *gIm, double *GRe, double *GIm)
{
std::vector<double> Gr2(m * n * 3);
std::vector<double> Gi2(m * n * 3); //temporary buffers
//calculate the fourier transform of the columns
for(int x = 0; x < n; x++)
for(int c = 0; c < 3; c++)
{
print(" % done",8, 0, RGB_White, 1);
print(50 * x / n, 0, 0, RGB_White, 1);
if(done()) end();
redraw();
//This is the 1D DFT:
for(int w = 0; w < m; w++)
{
Gr2[m * 3 * x + 3 * w + c] = Gi2[m * 3 * x + 3 * w + c] = 0;
for(int y = 0; y < m; y++)
{
double a= 2 * pi * w * y / float(m);
if(!inverse)a = a;
double ca = cos(a);
double sa = sin(a);
Gr2[m * 3 * x + 3 * w + c] += gRe[m * 3 * x + 3 * y + c] * ca  gIm[m * 3 * x + 3 * y + c] * sa;
Gi2[m * 3 * x + 3 * w + c] += gRe[m * 3 * x + 3 * y + c] * sa + gIm[m * 3 * x + 3 * y + c] * ca;
}
}
}
//calculate the fourier transform of the rows
for(int y = 0; y < m; y++)
for(int c = 0; c < 3; c++)
{
print(" % done",8, 0, RGB_White, 1);
print(50 + 50 * y / m, 0, 0, RGB_White, 1);
if(done()) end();
redraw();
//This is the 1D DFT:
for(int w = 0; w < n; w++)
{
GRe[m * 3 * w + 3 * y + c] = GIm[m * 3 * w + 3 * y + c] = 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[m * 3 * w + 3 * y + c] += Gr2[m * 3 * x + 3 * y + c] * ca  Gi2[m * 3 * x + 3 * y + c] * sa;
GIm[m * 3 * w + 3 * y + c] += Gr2[m * 3 * x + 3 * y + c] * sa + Gi2[m * 3 * x + 3 * y + c] * ca;
}
if(inverse)
{
GRe[m * 3 * w + 3 * y + c] /= n;
GIm[m * 3 * w + 3 * y + c] /= n;
}
else
{
GRe[m * 3 * w + 3 * y + c] /= m;
GIm[m * 3 * w + 3 * y + c] /= m;
}
}
}
}

void FFT2D(int n, int m, bool inverse, double *gRe, double *gIm, double *GRe, double *GIm)
{
int l2n = 0, p = 1; //l2n will become log_2(n)
while(p < n) {p *= 2; l2n++;}
int l2m = 0; p = 1; //l2m will become log_2(m)
while(p < m) {p *= 2; l2m++;}
m = 1 << l2m; n = 1 << l2n; //Make sure m and n will be powers of 2, otherwise you'll get in an infinite loop
//Erase all history of this array
for(int x = 0; x < m; x++) //for each column
for(int y = 0; y < m; y++) //for each row
for(int c = 0; c < 3; c++) //for each color component
{
GRe[3 * m * x + 3 * y + c] = gRe[3 * m * x + 3 * y + c];
GIm[3 * m * x + 3 * y + c] = gIm[3 * m * x + 3 * y + c];
}
//Bit reversal of each row
int j;
for(int y = 0; y < m; y++) //for each row
for(int c = 0; c < 3; c++) //for each color component
{
j = 0;
for(int i = 0; i < n  1; i++)
{
GRe[3 * m * i + 3 * y + c] = gRe[3 * m * j + 3 * y + c];
GIm[3 * m * i + 3 * y + c] = gIm[3 * m * j + 3 * y + c];
int k = n / 2;
while (k <= j) {j = k; k/= 2;}
j += k;
}
}
//Bit reversal of each column
double tx = 0, ty = 0;
for(int x = 0; x < n; x++) //for each column
for(int c = 0; c < 3; c++) //for each color component
{
j = 0;
for(int i = 0; i < m  1; i++)
{
if(i < j)
{
tx = GRe[3 * m * x + 3 * i + c];
ty = GIm[3 * m * x + 3 * i + c];
GRe[3 * m * x + 3 * i + c] = GRe[3 * m * x + 3 * j + c];
GIm[3 * m * x + 3 * i + c] = GIm[3 * m * x + 3 * j + c];
GRe[3 * m * x + 3 * j + c] = tx;
GIm[3 * m * x + 3 * j + c] = ty;
}
int k = m / 2;
while (k <= j) {j = k; k/= 2;}
j += k;
}
}
//Calculate the FFT of the columns
for(int x = 0; x < n; x++) //for each column
for(int c = 0; c < 3; c++) //for each color component
{
//This is the 1D FFT:
double ca = 1.0;
double sa = 0.0;
int l1 = 1, l2 = 1;
for(int l=0;l < l2n;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[3 * m * x + 3 * i1 + c]  u2 * GIm[3 * m * x + 3 * i1 + c];
double t2 = u1 * GIm[3 * m * x + 3 * i1 + c] + u2 * GRe[3 * m * x + 3 * i1 + c];
GRe[3 * m * x + 3 * i1 + c] = GRe[3 * m * x + 3 * i + c]  t1;
GIm[3 * m * x + 3 * i1 + c] = GIm[3 * m * x + 3 * i + c]  t2;
GRe[3 * m * x + 3 * i + c] += t1;
GIm[3 * m * x + 3 * i + c] += 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);
}
}
//Calculate the FFT of the rows
for(int y = 0; y < m; y++) //for each row
for(int c = 0; c < 3; c++) //for each color component
{
//This is the 1D FFT:
double ca = 1.0;
double sa = 0.0;
int l1= 1, l2 = 1;
for(int l = 0; l < l2m; 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[3 * m * i1 + 3 * y + c]  u2 * GIm[3 * m * i1 + 3 * y + c];
double t2 = u1 * GIm[3 * m * i1 + 3 * y + c] + u2 * GRe[3 * m * i1 + 3 * y + c];
GRe[3 * m * i1 + 3 * y + c] = GRe[3 * m * i + 3 * y + c]  t1;
GIm[3 * m * i1 + 3 * y + c] = GIm[3 * m * i + 3 * y + c]  t2;
GRe[3 * m * i + 3 * y + c] += t1;
GIm[3 * m * i + 3 * y + c] += 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);
}
}
int d;
if(inverse) d = n; else d = m;
for(int x = 0; x < n; x++) for(int y = 0; y < m; y++) for(int c = 0; c < 3; c++) //for every value of the buffers
{
GRe[3 * m * x + 3 * y + c] /= d;
GIm[3 * m * x + 3 * y + c] /= d;
}
}

//Draws an image
void draw(int xPos, int yPos, int n, int m, double *g, bool shift, bool neg128) //g is the image to be drawn
{
for(int x = 0; x < n; x++)
for(int y = 0; y < m; y++)
{
int x2 = x, y2 = y;
if(shift) {x2 = (x + n / 2) % n; y2 = (y + m / 2) % m;} //Shift corners to center
ColorRGB color;
//calculate color values out of the floating point buffer
color.r = int(g[3 * m * x2 + 3 * y2 + 0]);
color.g = int(g[3 * m * x2 + 3 * y2 + 1]);
color.b = int(g[3 * m * x2 + 3 * y2 + 2]);
if(neg128) color = color + RGB_Gray;
//negative colors give confusing effects so set them to 0
if(color.r < 0) color.r = 0;
if(color.g < 0) color.g = 0;
if(color.b < 0) color.b = 0;
//set color components higher than 255 to 255
if(color.r > 255) color.r = 255;
if(color.g > 255) color.g = 255;
if(color.b > 255) color.b = 255;
//plot the pixel
pset(x + xPos, y + yPos, color);
}
}

//Calculates the amplitude of *gRe and *gIm and puts the result in *gAmp
void calculateAmp(int n, int m, double *gAmp, double *gRe, double *gIm)
{
for(int x = 0; x < n; x++)
for(int y = 0; y < m; y++)
for(int c = 0; c < 3; c++)
{
gAmp[m * 3 * x + 3 * y + c] = sqrt(gRe[m * 3 * x + 3 * y + c] * gRe[m * 3 * x + 3 * y + c] + gIm[m * 3 * x + 3 * y + c] * gIm[m * 3 * x + 3 * y + c]);
}
}

//yeah, we're working with fixed sizes again...
const int N = 128; //the width of the image
const int M = 128; //the height of the image
double fRe[N][M][3], fIm[N][M][3], fAmp[N][M][3]; //the signal's real part, imaginary part, and amplitude
double FRe[N][M][3], FIm[N][M][3], FAmp[N][M][3]; //the FT's real part, imaginary part and amplitude
double fRe2[N][M][3], fIm2[N][M][3], fAmp2[N][M][3]; //will become the signal again after IDFT of the spectrum
double FRe2[N][M][3], FIm2[N][M][3], FAmp2[N][M][3]; //filtered spectrum
double pi = 3.1415926535897932384626433832795;
void draw(int xpos, int yPos, int n, int m, double *g, bool shift, bool neg128);
void FFT2D(int n, int m, bool inverse, double *gRe, double *gIm, double *GRe, double *GIm);
void calculateAmp(int n, int m, double *ga, double *gRe, double *gIm);

int main(int /*argc*/, char */*argv*/[])
{
screen(3 * N,4 * M + 8, 0,"2D FFT and Filters");
std::vector

if(keyPressed(SDLK_a)  changed)
{
for(int x = 0; x < N; x++)
for(int y = 0; y < M; y++)
for(int c = 0; c < 3; c++)
{
fRe2[x][y][c] = fRe[x][y][c];
} changed = 1;
} //no effect
if(keyPressed(SDLK_b))
{
for(int x = 0; x < N; x++)
for(int y = 0; y < M; y++)
for(int c = 0; c < 3; c++)
{
fRe2[x][y][c] = 255  fRe[x][y][c];
}
changed = 1;} //negative
if(keyPressed(SDLK_c))
{
for(int x = 0; x < N; x++)
for(int y = 0; y < M; y++)
for(int c = 0; c < 3; c++)
{
fRe2[x][y][c] = fRe[x][y][c]/2;
}
changed = 1;
} //half amplitude
//ETCETERA... *snip*
if(keyPressed(SDLK_SPACE)) endloop = 1;

if(changed)
{
//Draw the image and its FT
calculateAmp(N, M, fAmp2[0][0], fRe2[0][0], fIm2[0][0]);
draw(0, ytrans, N, M, fRe2[0][0], 0, 0); draw(N, 0+ytrans, N, M, fIm2[0][0], 0, 0); draw(2 * N, 0+ytrans, N, M, fAmp2[0][0], 0, 0); //draw real, imag and amplitude
FFT2D(N, M, 0, fRe2[0][0], fIm2[0][0], FRe[0][0], FIm[0][0]);
calculateAmp(N, M, FAmp[0][0], FRe[0][0], FIm[0][0]);
draw(0,M + ytrans, N, M, FRe[0][0], 1, 1); draw(N, M + ytrans, N, M, FIm[0][0], 1, 1); draw(2 * N, M + ytrans, N, M, FAmp[0][0], 1, 0); //draw real, imag and amplitude
print("Press az for effects and space to accept", 0, 0, RGB_White, 1, ColorRGB(128, 0, 0));
redraw();
}
changed = 0;
}
changed = 1;
while(!done())
{
readKeys();

if(keyPressed(SDLK_f))
{
for(int x = 0; x < N; x++)
for(int y = 0; y < M; y++)
for(int c = 0; c < 3; c++)
{
if(x * x + y * y < 64  (N  x) * (N  x) + (M  y) * (M  y) < 64  x * x + (M  y) * (M  y) < 64  (N  x) * (N  x) + y * y < 64)
{
FRe2[x][y][c] = FRe[x][y][c];
FIm2[x][y][c] = FIm[x][y][c];
}
else
{
FRe2[x][y][c] = 0;
FIm2[x][y][c] = 0;
}
}
changed = 1;
} //LP filter
if(keyPressed(SDLK_i))
{
for(int x = 0; x < N; x++)
for(int y = 0; y < M; y++)
for(int c = 0; c < 3; c++)
{
if(x * x + y * y > 16 && (N  x) * (N  x) + (M  y) * (M  y) > 16 && x * x + (M  y) * (M  y) > 16 && (N  x) * (N  x) + y * y > 16)
{
FRe2[x][y][c] = FRe[x][y][c];
FIm2[x][y][c] = FIm[x][y][c];
}
else
{
FRe2[x][y][c] = 0;
FIm2[x][y][c] = 0;
}
}
changed = 1;
} //HP filter
if(keyPressed(SDLK_s))
{
for(int x = 0; x < N; x++)
for(int y = 0; y < M; y++)
for(int c = 0; c < 3; c++)
{
if((x * x + y * y < 256  (N  x) * (N  x) + (M  y) * (M  y) < 256  x * x + (M  y) * (M  y) < 256  (N  x) * (N  x) + y * y < 256) && (x * x + y * y > 128 && (N  x) * (N  x) + (M  y) * (M  y) > 128 && x * x + (M  y) * (M  y) > 128 && (N  x) * (N  x) + y * y > 128))
{
FRe2[x][y][c] = FRe[x][y][c];
FIm2[x][y][c] = FIm[x][y][c];
}
else
{
FRe2[x][y][c] = 0;
FIm2[x][y][c] = 0;
}
}
changed = 1;
} //BP filter
//ETCETERA... *snip*

if(changed)
{
//after pressing a key: calculate the inverse!
FFT2D(N, M, 1, FRe2[0][0], FIm2[0][0], fRe[0][0], fIm[0][0]);
calculateAmp(N, M, fAmp[0][0], fRe[0][0], fIm[0][0]);
draw(0, 3 * M + ytrans, N, M, fRe[0][0], 0, 0); draw(N, 3 * M + ytrans, N, M, fIm[0][0], 0, 0); draw(2 * N, 3 * M + ytrans, N, M, fAmp[0][0], 0, 0);
calculateAmp(N, M, FAmp2[0][0], FRe2[0][0], FIm2[0][0]);
draw(0, 2 * M + ytrans, N, M, FRe2[0][0], 1, 1); draw(N, 2 * M + ytrans, N, M, FIm2[0][0], 1, 1); draw(2 * N, 2 * M + ytrans, N, M, FAmp2[0][0], 1, 0);
print("Press az for filters and arrows to change DC", 0, 0, RGB_White, 1, ColorRGB(128, 0, 0));
redraw();
}
changed=0;
}
return 0;
}
