Calculating the DFT in C++

When you learn about the Fourier transform and what it can show you about a signal, you immediately start thinking about its possible applications. The Fourier transform, however, deals with continuous time signals while, in practice, computers deal with discrete time signals (i.e. a sampled version of the original continuous time signal). When it comes to discrete time signal, you can calculate a discrete Fourier transform to get the frequency content of the signal.

Frequency and time domain representation of a 1 Hz sampled sine wave (sampling frequency is 8 Hz).

Typically, one should use a well known and tested library to perform this calculation for a number of reasons including accuracy and speed. However, it can be useful, for the sake of refreshing or improving one’s coding skills to try and code from scratch the algorithm for calculating the discrete Fourier transform. I decided to go this route to refresh and improve my C++ basic coding skills.

The DFT calculation is often implemented (for the sake of speed and efficiency) through the radix-2 algorithm or other forms of such algorithm. This allows to bring down the computation time to O(N log N) which makes a significant difference with respect to using the definition of the DFT, as N grows larger.

In my simple example, I will implement the definition of the DFT instead, that is, given a sequence of N samples $X$ the following formula is the definition of its DFT:

$$X_{k} = \sum_{n=0}^{N-1}x_{n}e^{-j2\pi k n/N}$$

Along with that, I implemented a lot of auxiliary functions used for printing vectors, saving data and output log information.

Let’s dive into the specifics of the program. The program is made of 2 files dft_main.cpp and dft_module.hpp: dft_main.cpp contains a set of variables that can be edited in order to obtain the desired output.

1. Input signal

As for the input, you can choose to use either a dummy signal (a 1 Hz sinewave) or an external signal which should be contained in a text file and be in the following format:

Header
sample1
sample2
…

Note that you can also have no header, provided that the variable “INPUT_FILE_HAS_HEADER” is set to false. By default it is assumed that the input file has a header. You can choose which input signal to use by changing the value of the variables “USE_EXTERNAL_DATA” and “INPUT_FILE_PATH”.

2. Output signal

If you want the output to be saved to a file, you should provide a path in “OUTPUT_FILE_PATH” and set “SAVE_TO_FILE” equal to true, otherwise the output will be simply printed to the console.

3. Sampling frequency

The sampling frequency can be set in Hz by changing the variable “FS”.

To run the program, simply compile it and then run it. You can use an IDE such as Code::Blocks. Finally, here is the code:

	/*
	===============================================================================
	Include statements
	===============================================================================
	*/
	#include<complex>
	#include<vector>
	#include<string>
	#include<dft_module.hpp>
	/*
	===============================================================================
	Namespaces statements
	===============================================================================
	*/
	using std::complex;
	using std::vector;
	using std::string;
	using dft_functions::calculate_dft;
	using dft_functions::get_number_of_lines;
	using dft_functions::print_info;
	using dft_functions::read_data_from_file;
	using dft_functions::save_dft_to_file;
	using dft_functions::test_dft_calculation;
	/*
	===============================================================================
	Constants declaration
	===============================================================================
	*/
	// External signal constants
	const double FS = 8.0; // Sampling frequency of the signal [Hz]
	const string OUTPUT_FILE_PATH = "C:/users/m/desktop/dft_calcs/data/test.txt"; // Output file path
	const string INPUT_FILE_PATH = "C:/users/m/desktop/dft_calcs/data/data.txt"; // Input file path
	const bool INPUT_FILE_HAS_HEADER = true; // Input file has header?
	const bool USE_EXTERNAL_DATA = false; // Use external data?
	// Save data to file?
	const bool SAVE_TO_FILE = false; // Save data to file?
	/*
	===============================================================================
	Main
	===============================================================================
	*/
	int main()
	{
	// Choose if use external data or sample data
	if(USE_EXTERNAL_DATA && (get_number_of_lines(INPUT_FILE_PATH, INPUT_FILE_HAS_HEADER) > 0) )
	{
	vector< complex<double> > signal = read_data_from_file(INPUT_FILE_PATH, INPUT_FILE_HAS_HEADER); // Sampled signal array
	int number_of_samples = signal.size();
	double frequency_resolution = FS / double(number_of_samples); // Frequency resolution of the DFT [Hz]
	vector<double> freq_axis(number_of_samples); // Frequency axis. From bins to frequency [Hz]
	for(int i = 0; i < number_of_samples; i++){ freq_axis[i] = double(i) * frequency_resolution; } // Fill the frequency axis
	// Compute DFT
	vector< complex<double> > dft = calculate_dft(signal);
	// Print info on the process and results
	print_info(number_of_samples, frequency_resolution, FS, USE_EXTERNAL_DATA, dft, freq_axis, INPUT_FILE_PATH);
	// Save to file
	if(SAVE_TO_FILE){ save_dft_to_file(dft, freq_axis, OUTPUT_FILE_PATH); }
	}else
	{
	// Run a dummy dft test calculation on a sinewave
	test_dft_calculation();
	}
	return 0;
	}

view raw dft_main.cpp hosted with ❤ by GitHub

	#include<iostream>
	#include<vector>
	#include<fstream>
	#include<complex>
	using std::cout;
	using std::cin;
	using std::endl;
	using std::vector;
	using std::complex;
	using std::string;
	using std::ofstream;
	using std::ifstream;
	namespace dft_constants
	{
	// General constants
	const complex<double> PI {{3.14159, 0}}; // PI
	const complex<double> IMAG_UNIT {{0, 1}}; // Imaginary unit
	// Test signal constants
	const double FS_TEST = 8.0; // Sampling frequency [Hz]
	const double TEST_FREQUENCY = 1.0; // Test signal frequency [Hz]
	const int NUMBER_OF_TEST_SAMPLES = 8; // Number of samples of the test signal (integer)
	}
	namespace dft_functions
	{
	// This function prints an array
	void print_array(const vector <complex<double> > &dft_array,
	const vector<double> &frequency_axis)
	{
	cout << endl << "Array of DFT: " << endl;
	cout << "F [Hz]" << "\t" << "Value" << endl;
	int number_of_samples = dft_array.size();
	for(int i=0; i < number_of_samples; i++)
	{
	cout << frequency_axis[i] << "\t" << dft_array[i] << endl;
	}
	}
	// This function prints info on the process
	void print_info(const int number_of_samples,
	const double &frequency_resolution,
	double fs,
	bool use_external_data,
	const vector <complex<double> > &dft_array,
	const vector<double> &frequency_axis,
	string input_file_path="")
	{
	if(use_external_data){ cout << "Computing DFT of data available from file " << input_file_path << endl; }
	else{ cout << "Computing DFT sample data: sine wave at 1 Hz" << endl; }
	cout << "Number of samples: " << number_of_samples << endl;
	cout << "Sampling frequency [Hz]: " << fs << endl;
	cout << "Frequency resolution [Hz]: " << frequency_resolution << endl;
	cout << "Maximum detectable frequency [Hz]: " << fs/2.0 << endl << endl;
	print_array(dft_array, frequency_axis);
	}
	// This function calculates the normalized DFT
	vector< complex<double> > calculate_dft(const vector< complex<double> > &signal)
	{
	vector< complex<double> > dft(signal.size()); // DFT vector
	double N = double(signal.size()); // Number of samples
	complex<double> temp = {{0, 0}}; // Temporary loop variable
	for(int k = 0; k < signal.size(); k++)
	{
	for(int n = 0; n < signal.size(); n++)
	{
	temp = {{double(-1*2*k*n) / N, 0}};
	dft[k] += signal[n] * exp(dft_constants::IMAG_UNIT * dft_constants::PI * temp) / N;
	}
	}
	return dft;
	}
	// This function checks if a file exists
	bool check_if_file_exists(string file_path)
	{
	ifstream file(file_path);
	bool file_existance_status = file.good();
	file.close();
	return file_existance_status;
	}
	// This function saves a complex array to a file
	void save_dft_to_file(const vector< complex<double> > &dft,
	const vector<double> &freq_axis,
	string file_path)
	{
	// Check if the file already exists
	if(check_if_file_exists(file_path))
	{
	cout << "File: " << file_path << " Already exists... OVERWRITE? (y/n) ";
	string answer;
	cin >> answer;
	if(answer != "y")
	{
	cout << "Data not saved!" << endl;
	return;
	}
	cout << "Overwriting..." << endl;
	}
	// Open file
	ofstream file(file_path);
	// If the file is open save data to it
	if(file.is_open())
	{
	// Set header
	file << "Frequency" << "," << "real_part" << "," << "imaginary_part" << "\n";
	// Feed data to the file
	for(int i=0; i < dft.size(); i++)
	{
	file << freq_axis[i] << "," << dft[i].real() << "," << dft[i].imag() << "\n";
	}
	// Close stream
	file.close();
	cout << "Data successfully saved to file in the following path: " << file_path << endl;
	}else
	{
	cout << "Unable to open file... Data has not been saved." << endl;
	}
	}
	// This function gets the number of lines in a file (less the header line)
	int get_number_of_lines(string file_path, bool header)
	{
	// Open file
	ifstream file(file_path);
	int number_of_lines = 0;
	string line;
	// Count the number of lines (less the header)
	if(file.is_open())
	{
	while(getline(file, line))
	{
	number_of_lines ++;
	}
	if(header)
	{
	number_of_lines --;
	}
	}
	// Close file
	file.close();
	return number_of_lines;
	}
	// This function reads and parses a single column csv file
	vector< complex<double> > read_data_from_file(string file_path, bool header)
	{
	// Get the number of lines
	int number_of_lines = get_number_of_lines(file_path, header);
	// Open the file stream
	ifstream file(file_path);
	// Read file and load data into a vector
	if(file.is_open() && (number_of_lines > 0))
	{
	vector< complex<double> > output_vector(number_of_lines);
	string line = "";
	string splitted_string = "";
	string delimiter = "\n";
	int i=0;
	if(header)
	{
	getline(file, line);
	cout << "Header removed!" << endl;
	}
	while( getline(file, line) )
	{
	splitted_string = line.substr(0, line.find(delimiter));
	// Store the splitted string as a double using stod()
	output_vector[i] = {stod(splitted_string), 0};
	i++;
	}
	file.close();
	cout << "File opened and parsed successfully!" << endl;
	return output_vector;
	}else
	{
	cout << "Cannot open file!" << endl;
	vector< complex<double> > output_vector {{0, 0}};
	return output_vector;
	}
	}
	// Run a dummy DFT calculation on a sinewave of TEST_FREQUENCY Hz sampled at FS_TEST Hz.
	void test_dft_calculation()
	{
	vector< complex<double> > signal (dft_constants::NUMBER_OF_TEST_SAMPLES); // Sampled signal array
	vector<double> freq_axis(dft_constants::NUMBER_OF_TEST_SAMPLES); // Frequency axis. From bins to frequency [Hz]
	double frequency_resolution = dft_constants::FS_TEST/double(dft_constants::NUMBER_OF_TEST_SAMPLES); // Frequency resolution of the DFT [Hz]
	// Build signal and frequency axis
	for(int i = 0; i < dft_constants::NUMBER_OF_TEST_SAMPLES; i++)
	{
	freq_axis[i] = double(i) * frequency_resolution;
	signal[i] = {sin(2 * dft_constants::PI.real() * dft_constants::TEST_FREQUENCY * (1/dft_constants::FS_TEST) * double(i)), 0};
	}
	// Compute DFT
	vector< complex<double> > dft = calculate_dft(signal);
	// Print info on the process and results
	print_info(dft_constants::NUMBER_OF_TEST_SAMPLES, frequency_resolution, dft_constants::FS_TEST, false, dft, freq_axis);
	}
	}

view raw dft_module.hpp hosted with ❤ by GitHub