Faculty of Engineering: C++ Program for DFT & Signal Processing

This assignment showcases several C++ programs designed to perform a Discrete Fourier Transform (DFT) on input signals. The programs read signal values, calculate the DFT using different approaches, including brute force and twiddle factors, and output the results to a MATLAB/Octave-compatible file for plotting. The DFT process involves calculating real and imaginary parts of the DFT coefficients to identify frequency components. The code includes signal generation, DFT calculation, and power spectrum computation, with the final output visualized through plots of the original signal, real and imaginary parts of the DFT, and the power spectrum. Desklib provides this assignment as a resource for students studying signal processing and DFT implementations.

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;

}