Program for DFT process

This program calculates the Discrete Fourier Transform (DFT) of a given signal using twiddle factors and brute force methods in C++. It generates a random discrete-time signal and calculates the DFT and power spectrum up to the Nyquist frequency. The output is saved in a MATLAB/Octave M-file for plotting.

Program for DFT process
(1)
#include <iostream .h>
#include <iomanip.h>
#include <math.h>
// N is assumed to be greater than 4 and a power of 2
#define N 360
#define PI2 6.2832
// Twiddle factors (360th roots of unity)
const float W[] = {
0.0175, 0.0349, 0.0523, 0.0698, 0.0872, 0.1045, 0.1219, 0.1392, 0.1564, 0.1736,
0.1908, 0.2079, 0.2250, 0.2419, 0.2588, 0.2756, 0.2924, 0.3090, 0.3256, 0.3420,
0.3584, 0.3746, 0.3907, 0.4067, 0.4226, 0.4384, 0.4540, 0.469, 0.4848, 0.5000,
0.5150, 0.5299, 0.5446, 0.5592, 0.5736, 0.5878, 0.6018, 0.6157 0.6293, 0.6428,
0.6561, 0.6691, 0.6820, 0.6947, 0.7071, 0.7193, 0.731, 0.7431, 0.7547, 0.7660,
0.7771 0.7880, 0.7986, 0.8090, 0.8192, 0.8290, 0.838, 0.8480, 0.8572, 0.8660,
0.8746, 0.8829, 0.8910, 0.8988, 0.9063, 0.9135, 0.9205, 0.9272, 0.9336, 0.9397,
0.9455, 0.9511, 0.9563, 0.9613, 0.9659, 0.9703, 0.9744, 0.9781, 0.9816, 0.9848,
0.9877, 0.9903, 0.9925, 0.9945, 0.9962, 0.9976, 0.9986, 0.9994, 0.9998, 1.0000,
0.9998, 0.9994, 0.9986, 0.9976, 0.9962, 0.9945, 0.9925, 0.9903, 0.9877, 0.9848,
0.9816, 0.9781, 0.9744, 0.9703, 0.9659, 0.9613, 0.9563, 0.9511, 0.9455, 0.9397,
0.9336, 0.9272, 0.9205, 0.9135, 0.9063, 0.8988, 0.8910, 0.8829, 0.8746, 0.8660,
0.8572, 0.8480, 0.8387, 0.8290, 0.8192, 0.8090, 0.7986, 0.7880, 0.7771, 0.7660,
0.7547, 0.7431, 0.7314, 0.7193, 0.7071, 0.6947, 0.6691, 0.6820, 0.561, 0.6428,
0.6293, 0.6157, 0.6018, 0.5878, 0.5736, 0.5592, 0.5446, 0.5299, 0.5150, 0.5000,
0.4848, 0.4695, 0.4540, 0.4384, 0.4226, 0.40677, 0.3907, 0.3746, 0.3584, 0.3420,
0.3256, 0.3090, 0.2924, 0.2588, 0.2756, 0.2419, 0.2250, 0.2079, 0.1908, 0.1736,
0.1564, 0.1392, 0.1219, 0.1045, 0.0698, 0.0523, 0.0349, 0.0175, 0.0000,-0.0175,
-0.0349,-0.0523,-0.0698,-0.0872, 0.104, -0.1219,-0.1392, -0.1564,-0.1736,-0.1908,
-0.2079,-0.2250,-0.2419,-0.2756,-0.3090,-0.3256, 0.3584,-0.3746, -0.3907, -0.4067,
-0.4226,-0.4384,-0.4540, -0.4695,-0.4848,-0.5000,-0.5150, -0.5299,-0.5446,-0.5592,

-0.5736,-0.5878,-0.6018,-0.6157,-0.6293,-0.6428,-0.6561,-0.6691 -0.6820,-0.6947,
0.7071,-0.7193,-0.7314,-0.7431, -0.7547,-0.7660,-0.7771, -0.7880, -0.7986, 0.8090,
-0.8192,-0.8290,-0.8387,-0.8480,-0.8572,-0.8660,-0.8746,-0.8829,-0.8910,-0.8988,
-0.9063,-0.9135, -0.9205,-0.9272,-0.9336,-0.9397,-0.9455,-0.9511,-0.9563,-0.9613,
-0.9659,-0.9703,-0.9744,-0.9781, 0.9816,-0.9848,-0.9877,-0.9903,-0.9925,-0.9945,
-0.9962,-0.9976,-0.9986,-0.9994, 0.9998, 1.0000,-0.9998,-0.9994,-0.9986,-0.9976,
-0.9962,-0.9945,-0.9925,-0.9903,-0.9877,-0.9848, 0.9816,-0.9781,-0.9744,-0.9703,
-0.9659,-0.9613,-0.9563,-0.9511,-0.9455,-0.9397,-0.9336, 0.9272,-0.9205,-0.9135,
-0.9063,-0.8988,-0.8910,-0.8829,-0.8746,-0.8660,-0.8572,-0.8480, 0.8387,-0.8290,
-0.8192,-0.8090,-0.7986,-0.7880,-0.7771,-0.7660,-0.7547,-0.7431,-0.7314,-0.7193,
-0.7071,-0.6947,-0.6820,-0.6691,-0.6561,-0.6428,-0.6293,-0.6157,-0.6018,-0.5878,
-0.5736, 0.5592,-0.5446,-0.5299,-0.5150,-0.5000,-0.4848,-0.4695,-0.4540,-0.4384,
-0.4226,-0.4067,-0.3907,-0.3746,-0.3584, -0.3420, -0.3256, -0.3090, -0.2924,-0.2756,
-0.2588, -0.2419, -0.2250, -0.2079, -0.1908, -0.1736, -0.1564, -0.1392, -0.1219, -0.1045,
-0.0872, -0.0698,-0.0523, -0.0349, 0.0175, 0.0000
};
int main()
{
// time and frequency domain data arrays
int n, k; // time and frequency domain indices
float x[N]; // discrete-time signal, x
float Xre[N/2+1], Xim[N/2+1]; // DFT of x (real and imaginary parts)
float P[N/2+1]; // power spectrum of x
// Generate random discrete-time signal x in range (-1,+1)
srand(time(0));
for (n=0 ; n<N ; ++n) x[n] = ((2.0 * rand()) / RAND_MAX) - 1.0 + sin(PI2 * n * 5.7 / N);
// Calculate DFT and power spectrum up to Nyquist frequency
int to_sin = 3*N/4; // index offset for sin
int a, b;
for (k=0 ; k<=N/2 ; ++k)
{

Xre[k] = 0; Xim[k] = 0;
a = 0; b = to_sin;
for (n=0 ; n<N ; ++n)
{
Xre[k] += x[n] * W[a%N];
Xim[k] -= x[n] * W[b%N];
a += k; b += k;
}
P[k] = Xre[k]*Xre[k] + Xim[k]*Xim[k];
}
// Output results to MATLAB / Octave M-file for plotting
FILE *f = fopen("dftplots.m", "w");
fprintf(f, "n = [0:%d];\n", N-1);
fprintf(f, "x = [ ");
for (n=0 ; n<N ; ++n) fprintf(f, "%f ", x[n]);
fprintf(f, "];\n");
fprintf(f, "Xre = [ ");
for (k=0 ; k<=N/2 ; ++k) fprintf(f, "%f ", Xre[k]);
fprintf(f, "];\n");
fprintf(f, "Xim = [ ");
for (k=0 ; k<=N/2 ; ++k) fprintf(f, "%f ", Xim[k]);
fprintf(f, "];\n");
fprintf(f, "P = [ ");
for (k=0 ; k<=N/2 ; ++k) fprintf(f, "%f ", P[k]);
fprintf(f, "];\n");
fprintf(f, "subplot(3,1,1)\nplot(n,x)\n");
fprintf(f, "xlim([0 %d])\n", N-1);
fprintf(f, "subplot(3,1,2)\nplot([0:%d],Xre,[0:%d],Xim)\n", N/2, N/2);
fprintf(f, "xlim([0 %d])\n", N/2);
fprintf(f, "subplot(3,1,3)\nstem([0:%d],P)\n", N/2);
fprintf(f, "xlim([0 %d])\n", N/2);
fclose(f);

// exit normally
return 0;}
(2)
#include <iostream.h>
#include <iomanip.h>
#include <math.h>
#define n=0:-0.5000
#define N =360
#define PI2 6.2832
int main()
{
// time and frequency domain data arrays
int n, k; // indices for time and frequency domains
float x[N]; // discrete-time signal, x
float Xre[N], Xim[N]; // DFT of x (real and imaginary parts)
float P[N]; // power spectrum of x
// Generate random discrete-time signal x in range (-1,+1)
srand(time(0));
for (n=0 ; n<N ; ++n) x[n] = ((2.0 * rand()) / RAND_MAX) - 1.0;
// Calculate DFT of x using brute force
for (k=0 ; k<N ; ++k)
{
// Real part of X[k]
Xre[k] = 0;
for (n=0 ; n<N ; ++n) Xre[k] += x[n] * cos(n * k * PI2 / N);
// Imaginary part of X[k]
Xim[k] = 0;
for (n=0 ; n<N ; ++n) Xim[k] -= x[n] * sin(n * k * PI2 / N);

// Power at kth frequency bin
P[k] = Xre[k]*Xre[k] + Xim[k]*Xim[k];
}
// Output results to MATLAB / Octave M-file for plotting
FILE *f = fopen("dftplots.m", "w");
fprintf(f, "n = [0:%d];\n", N-1);
fprintf(f, "x = [ ");
for (n=0 ; n<N ; ++n) fprintf(f, "%f ", x[n]);
fprintf(f, "];\n");
fprintf(f, "Xre = [ ");
for (k=0 ; k<N ; ++k) fprintf(f, "%f ", Xre[k]);
fprintf(f, "];\n");
fprintf(f, "Xim = [ ");
for (k=0 ; k<N ; ++k) fprintf(f, "%f ", Xim[k]);
fprintf(f, "];\n");
fprintf(f, "P = [ ");
for (k=0 ; k<N ; ++k) fprintf(f, "%f ", P[k]);
fprintf(f, "];\n");
fprintf(f, "subplot(3,1,1)\nplot(n,x)\n");
fprintf(f, "xlim([0 %d])\n", N-1);
fprintf(f, "subplot(3,1,2)\nplot(n,Xre,n,Xim)\n");
fprintf(f, "xlim([0 %d])\n", N-1);
fprintf(f, "subplot(3,1,3)\nstem(n,P)\n");
fprintf(f, "xlim([0 %d])\n", N-1);
fclose(f);
// exit normally
return 0;
}

(3)
#include <stdio.h>
#include <iostream .h>
#include <math.h>
// N is assumed to be greater than 4 and a power of 2
#define N =360
#define PI2 6.2832
// Twiddle factors (360th roots of unity)
const float W[] = {
-0.8177, -0.8010,
-0.7832, -0.7644, -0.7446, -0.7238, -0.7019, -0.6791, -0.6554, -0.6307, -0.6051, -0.5787,
-0.5514, -0.5234, -0.4945, -0.4650, -0.4347, -0.4038, -0.3724, -0.3403, -0.3078, -0.2748,
-0.2413, -0.2075, -0.1734, -0.1391, -0.1045, -0.0697, -0.0349, 0.0000, 0.0349, 0.0697,
0.1045, 0.1391, 0.1734, 0.2075, 0.2413, 0.2748, 0.3078, 0.3403, 0.3724, 0.4038,
0.4347, 0.4650, 0.4945, 0.5234, 0.5514, 0.5787, 0.6051, 0.6307, 0.6554, 0.6791,
0.7019, 0.7238, 0.7446, 0.7644, 0.7832, 0.8010, 0.8177, 0.8333, 0.8479, 0.8614,
0.8739, 0.8852, 0.8956, 0.9048, 0.9130, 0.9202, 0.9263, 0.9315, 0.9356, 0.9388,
0.9410, 0.9424, 0.9428, 0.9424, 0.9411, 0.9391, 0.9363, 0.9327, 0.9285, 0.9236,
0.9181, 0.9120, 0.9054, 0.8983, 0.8908, 0.8829, 0.8746, 0.8660, 0.8572, 0.8481,
0.8389, 0.8295, 0.8200, 0.8105, 0.8010, 0.7916, 0.7823, 0.7730, 0.7640, 0.7551,
0.7465, 0.7382, 0.7302, 0.7226, 0.7153, 0.7085, 0.7021, 0.6961, 0.6907, 0.6858,
0.6814, 0.6775, 0.6742, 0.6715, 0.6694, 0.6679, 0.6670, 0.6667, 0.6670, 0.6679,
0.6694, 0.6715, 0.6742, 0.6775, 0.6814, 0.6858, 0.6907, 0.6961, 0.7021, 0.7085,
0.7153, 0.7226, 0.7302, 0.7382, 0.7465, 0.7551, 0.7640, 0.7730, 0.7823, 0.7916,
0.8010, 0.8105, 0.8200, 0.8389, 0.8481, 0.8572, 0.8660, 0.8746, 0.8829, 0.8908,
0.8983, 0.9054, 0.9120, 0.9181, 0.9236, 0.9285, 0.9327, 0.9363, 0.9391, 0.9411,
0.9424, 0.9428, 0.9424, 0.9410, 0.9388, 0.9356, 0.9315, 0.9263, 0.9202, 0.9130,
0.9048, 0.8956, 0.8852, 0.8739, 0.8614, 0.8479, 0.8333, 0.8177, 0.8010, 0.7832,
0.7644, 0.7446, 0.7238, 0.7019, 0.6791, 0.6554, 0.6307, 0.6051, 0.5787, 0.5514,
0.5234, 0.4945, 0.4650, 0.4347, 0.4038, 0.3724, 0.3403, 0.3078, 0.2748, 0.2413,
0.2075, 0.1734, 0.1391, 0.1045, 0.0697, 0.0349, 0.0000, -0.0349,-0.0697, -0.1045,
-0.1391, -0.1734, -0.2075, -0.2413, -0.2748, -0.3078, -0.3403, -0.3724, -0.4038, -0.4347,

-0.4650, -0.4945, -0.5234, -0.5514, -0.5787, -0.6051, -0.6307, -0.6554, -0.6791, -0.7019,
-0.7238, -0.7446, -0.7644, -0.7832, -0.8010, -0.8177, -0.8333, -0.8479, -0.8614,-0.8739,
-0.9356, -0.9388, -0.9410, -0.9424, -0.9428, -0.9424, -0.9411, -0.9391, -0.9363, -0.9327,
-0.9285, -0.9236, -0.9181,-0.9120, -0.9054, -0.8983, -0.8908, -0.8829 -0.8746, -0.8660,
-0.8572, -0.8481, -0.8389, -0.8295, -0.8200, -0.8105, -0.8010, -0.7916, -0.7823, -0.7730,
-0.7640, -0.7551, -0.7465, -0.7382, -0.7302, -0.7226, -0.7153, -0.7085, -0.7021, -0.6961,
-0.6907, -0.6858, -0.6814, -0.6775, -0.6742, -0.6715, -0.6694, -0.6679, -0.6670, -0.6667,
-0.6670, -0.6679, -0.6694, -0.6715, -0.6742, -0.6775, -0.6814, -0.6858, -0.6907, -0.6961,
-0.7021, -0.7085, -0.7153, -0.7226, -0.7302, -0.7382, -0.7465, -0.7551, -0.7640, -0.7730,
-0.7823,-0.7916,-0.8010,-0.8105, -0.8200, -0.8295, -0.8389, -0.8481, -0.8572, -0.8660
-0.8746, -0.8829, -0.8908, -0.8983, -0.9054, -0.9120, -0.9181, -0.9236, -0.9285, -0.9327
-0.9363, -0.9391,-0.9411, -0.9424, -0.9428, -0.9424, -0.9410, -0.9388, -0.9356, -0.9315,
-0.9263, -0.9202, -0.9130, -0.9048, -0.8956, -0.8852, -0.8739, -0.8614, -0.8479, -0.8333
};
int main()
{
// time and frequency domain data arrays
int n, k; // time and frequency domain indices
float x[N]; // discrete-time signal, x
float Xre[N/2+1], Xim[N/2+1]; // DFT of x (real and imaginary parts)
float P[N/2+1]; // power spectrum of x
// Generate random discrete-time signal x in range (-1,+1)
srand(time(0));
for (n=0 ; n<N ; ++n) x[n] = ((2.0 * rand()) / RAND_MAX) - 1.0 + sin(PI2 * n * 5.7 / N);
// Calculate DFT and power spectrum up to Nyquist frequency
int to_sin = 3*N/4; // index offset for sin
int a, b;
for (k=0 ; k<=N/2 ; ++k)
{
Xre[k] = 0; Xim[k] = 0;

a = 0; b = to_sin;
for (n=0 ; n<N ; ++n)
{
Xre[k] += x[n] * W[a%N];
Xim[k] -= x[n] * W[b%N];
a += k; b += k;
}
P[k] = Xre[k]*Xre[k] + Xim[k]*Xim[k];
}
// Output results to MATLAB / Octave M-file for plotting
FILE *f = fopen("dftplots.m", "w");
fprintf(f, "n = [0:%d];\n", N-1);
fprintf(f, "x = [ ");
for (n=0 ; n<N ; ++n) fprintf(f, "%f ", x[n]);
fprintf(f, "];\n");
fprintf(f, "Xre = [ ");
for (k=0 ; k<=N/2 ; ++k) fprintf(f, "%f ", Xre[k]);
fprintf(f, "];\n");
fprintf(f, "Xim = [ ");
for (k=0 ; k<=N/2 ; ++k) fprintf(f, "%f ", Xim[k]);
fprintf(f, "];\n");
fprintf(f, "P = [ ");
for (k=0 ; k<=N/2 ; ++k) fprintf(f, "%f ", P[k]);
fprintf(f, "];\n");
fprintf(f, "subplot(3,1,1)\nplot(n,x)\n");
fprintf(f, "xlim([0 %d])\n", N-1);
fprintf(f, "subplot(3,1,2)\nplot([0:%d],Xre,[0:%d],Xim)\n", N/2, N/2);
fprintf(f, "xlim([0 %d])\n", N/2);
fprintf(f, "subplot(3,1,3)\nstem([0:%d],P)\n", N/2);
fprintf(f, "xlim([0 %d])\n", N/2);
fclose(f);
// exit normally

return 0;
}
(4)
#include <iostream.h>
#include <iomanip.h>
#include <math.h>
#define n=0:-0.8434
#define N =360
#define PI2 6.2832
int main()
{
// time and frequency domain data arrays
int n, k; // indices for time and frequency domains
float x[N]; // discrete-time signal, x
float Xre[N], Xim[N]; // DFT of x (real and imaginary parts)
float P[N]; // power spectrum of x
// Generate random discrete-time signal x in range (-1,+1)
srand(time(0));
for (n=0 ; n<N ; ++n) x[n] = ((2.0 * rand()) / RAND_MAX) - 1.0;
// Calculate DFT of x using brute force
for (k=0 ; k<N ; ++k)
{
// Real part of X[k]
Xre[k] = 0;
for (n=0 ; n<N ; ++n) Xre[k] += x[n] * cos(n * k * PI2 / N);

// Imaginary part of X[k]
Xim[k] = 0;
for (n=0 ; n<N ; ++n) Xim[k] -= x[n] * sin(n * k * PI2 / N);
// Power at kth frequency bin
P[k] = Xre[k]*Xre[k] + Xim[k]*Xim[k];
}
// Output results to MATLAB / Octave M-file for plotting
FILE *f = fopen("dftplots.m", "w");
fprintf(f, "n = [0:%d];\n", N-1);
fprintf(f, "x = [ ");
for (n=0 ; n<N ; ++n) fprintf(f, "%f ", x[n]);
fprintf(f, "];\n");
fprintf(f, "Xre = [ ");
for (k=0 ; k<N ; ++k) fprintf(f, "%f ", Xre[k]);
fprintf(f, "];\n");
fprintf(f, "Xim = [ ");
for (k=0 ; k<N ; ++k) fprintf(f, "%f ", Xim[k]);
fprintf(f, "];\n");
fprintf(f, "P = [ ");
for (k=0 ; k<N ; ++k) fprintf(f, "%f ", P[k]);
fprintf(f, "];\n");
fprintf(f, "subplot(3,1,1)\nplot(n,x)\n");
fprintf(f, "xlim([0 %d])\n", N-1);
fprintf(f, "subplot(3,1,2)\nplot(n,Xre,n,Xim)\n");
fprintf(f, "xlim([0 %d])\n", N-1);
fprintf(f, "subplot(3,1,3)\nstem(n,P)\n");
fprintf(f, "xlim([0 %d])\n", N-1);
fclose(f);
// exit normally
return 0;

}