Commit 7d9e1adf authored by Vyomkesh's avatar Vyomkesh
Browse files

Initial commit

parents
# Project exclude paths
/cmake-build-debug/
\ No newline at end of file
# Default ignored files
/shelf/
/workspace.xml
# Datasource local storage ignored files
/dataSources/
/dataSources.local.xml
# Editor-based HTTP Client requests
/httpRequests/
<?xml version="1.0" encoding="UTF-8"?>
<module classpath="CMake" type="CPP_MODULE" version="4" />
\ No newline at end of file
<?xml version="1.0" encoding="UTF-8"?>
<project version="4">
<component name="CMakeWorkspace" PROJECT_DIR="$PROJECT_DIR$" />
<component name="CidrRootsConfiguration">
<sourceRoots>
<file path="$PROJECT_DIR$/headers" />
<file path="$PROJECT_DIR$/sources" />
</sourceRoots>
</component>
</project>
\ No newline at end of file
<?xml version="1.0" encoding="UTF-8"?>
<project version="4">
<component name="ProjectModuleManager">
<modules>
<module fileurl="file://$PROJECT_DIR$/.idea/fast_fourier.iml" filepath="$PROJECT_DIR$/.idea/fast_fourier.iml" />
</modules>
</component>
</project>
\ No newline at end of file
cmake_minimum_required(VERSION 3.17)
project(fast_fourier)
set(CMAKE_CXX_STANDARD 14)
set(SOURCE_FILES sources/fast_fourier.cpp headers/fast_fourier.h)
add_library(fast_fourier_lib SHARED ${SOURCE_FILES})
find_package(OpenMP)
if (OpenMP_CXX_FOUND)
target_link_libraries(fast_fourier_lib OpenMP::OpenMP_CXX)
endif ()
add_executable(fast_fourier main.cpp )
TARGET_LINK_LIBRARIES(fast_fourier fast_fourier_lib)
//
// Created by vyomkesh jha on 27.12.2020.
//
#ifndef FAST_FOURIER_FAST_FOURIER_H
#define FAST_FOURIER_FAST_FOURIER_H
#include <stdlib.h>
#include <iostream>
#include <string.h>
#include <stdio.h>
#include <time.h>
#include <math.h>
#include <omp.h>
#include <complex>
using namespace std;
class fast_fourier {
public:
void perform_bit_reversal(int n, double *xReal, double *xImg, double log2_of_n);
void fill_array(int n, double *xReal, double *xImg, double *yReal, double *yImg);
void calculate_fft_par(int log2_of_n, double *xReal, double *xImg, double *wReal, double *wImg, int N);
timespec get_time_taken(timespec s, timespec f);
void get_twiddle_factors(int n, int N, double *wReal, double *wImg);
};
#endif //FAST_FOURIER_FAST_FOURIER_H
//
// Created by vyomkesh jha on 27.12.2020.
//
#include "headers/fast_fourier.h"
int main(int argc, char** argv)
{
fast_fourier *transformer = new fast_fourier();
if (argc != 3)
{
cerr << "usage: FFT <size of array> <number of threads>" << endl;
exit(1);
}
int n = atoi(argv[1]); //Assign size of array
if (!(n > 0) || !((n & (n - 1)) == 0)) //Test if
{
cerr << "usage: <size of array> must be power of 2" << endl;
exit(1);
}
int p = atoi(argv[2]); //Assign number of threads to use
double log2_of_n = log10(n) / log10(2); //Set for use later
double *xReal = new double[n];
double *xImg = new double[n];
double *yReal = new double[n];
double *yImg = new double[n];
double *wReal = new double[(int) log2_of_n];
double *wImg = new double[(int) log2_of_n];
double *W_OracleR = new double[n / 2];
double *W_OracleI = new double[n / 2];
omp_set_num_threads(p); //Set how many proc's to user
transformer->fill_array(n, xReal, xImg, yReal, yImg);
timespec sTime;
clock_gettime(CLOCK_MONOTONIC, &sTime); //Set start time
transformer->perform_bit_reversal(n, xReal, xImg, log2_of_n);
transformer->get_twiddle_factors(log2_of_n, n, wReal, wImg);
transformer->calculate_fft_par(log2_of_n, xReal, xImg, wReal, wImg, n);
timespec fTime;
clock_gettime(CLOCK_MONOTONIC, &fTime); //Set end time
double frac = 1.0 / (double) n; //Set fraction var
double error = 0.0;
for (int i = 0; i < n; i++)
{ //Iterate through values and compile error with FFT ( FFT ( X(1:N) ) ) == N * X(1:N)
//Thus we want close to zero as possible
error = error + pow(yReal[i] - frac * xReal[i], 2) + pow(yImg[i] - frac * xImg[i], 2);
}
error = sqrt(frac * error); //Calculate final error
cout << "Parallel FFT Time-elapsed: " << transformer->get_time_taken(sTime, fTime).tv_sec
<< "." << transformer->get_time_taken(sTime, fTime).tv_nsec << " seconds; ";// << endl;
cout << "Error (Mine): " << error << endl;
return 0;
}
//
// Created by vyomkesh jha on 27.12.2020.
//
#include "../headers/fast_fourier.h"
void fast_fourier::perform_bit_reversal(int n, double *xReal, double *xImg, double log2_of_n) {
int i;
double temp;
#pragma omp parallel for default(none) private(i, temp) shared(n, xReal, xImg, log2_of_n)
for (i = 1; i < n; i++) {
int k, rev = 0;
int inp = i;
for (k = 0; k < log2_of_n; k++) {
rev = (rev << 1) | (inp & 1);
inp >>= 1;
}
if (rev <= i) continue; //Skip if already done
temp = xReal[i]; //Store into temp values and swap
xReal[i] = xReal[rev];
xReal[rev] = temp;
temp = xImg[i];
xImg[i] = xImg[rev];
xImg[rev] = temp;
}
}
void fast_fourier::fill_array(int n, double *xReal, double *xImg, double *yReal, double *yImg) {
int i;
struct drand48_data buffer;
srand48_r(time(NULL) ^ omp_get_thread_num(), &buffer); //Seed for each thread
//#pragma omp parallel for default(none) private(i) shared(x,n)
for (i = 0; i < n; i++) {
drand48_r(&buffer, &xReal[i]); //Store random double into xReal
yReal[i] = xReal[i]; //Copy
drand48_r(&buffer, &xImg[i]); //Store random double into xImg
}
}
void fast_fourier::calculate_fft_par(int log2_of_n, double *xReal, double *xImg, double *wReal, double *wImg, int N) {
int n, d, i, k;
double temp_w, temp_x;
for (n = 1; n < log2_of_n + 1; n++) { //Iterate through time slice
d = pow(2, n);
#pragma omp parallel for default(none) private(k, i, temp_w, temp_x) shared(n, N, d, xReal, xImg, wReal, wImg, log2_of_n)
for (k = 0; k < (d / 2); k++) { //Iterate through even and odd elements
for (i = k; i < N; i += d) { //Butterfly operation
temp_w = wReal[n - 1] * xReal[i + (d / 2)]; //Multiply by twiddle factor
temp_x = xReal[i]; //Store in temp's and restore with correct values
xReal[i] = temp_w + temp_x;
xReal[i + (d / 2)] = temp_x - temp_w;
temp_w = wImg[n - 1] * xImg[i + (d / 2)]; //Multiply by twiddle factor
temp_x = xImg[i]; //Store in temp's and restore with correct values
xImg[i] = temp_w + temp_x;
xImg[i + (d / 2)] = temp_x - temp_w;
}
}
}
}
timespec fast_fourier::get_time_taken(timespec s, timespec f) {
timespec totTime;
if ((f.tv_nsec - s.tv_nsec) < 0) //so we do not return a negative incorrect value
{
totTime.tv_sec = f.tv_sec - s.tv_sec - 1;
totTime.tv_nsec = 1000000000 + f.tv_nsec - s.tv_nsec;
} else {
totTime.tv_sec = f.tv_sec - s.tv_sec;
totTime.tv_nsec = f.tv_nsec - s.tv_nsec;
}
return totTime;
}
void fast_fourier::get_twiddle_factors(int log2_of_n, int n, double wReal[], double wImg[]) {
int i;
#pragma omp parallel for default(none) private(i) shared(n, wReal, wImg, log2_of_n)
for (i = 0; i < n; i++) {
wReal[i] = cos(((double) i * 2.0 * M_PI) / ((double) n));
wImg[i] = sin(((double) i * 2.0 * M_PI) / ((double) n));
}
}
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment