2. Find or write (think about recursion, here) a program to calculate FFT. Comment the program in your own words.
#include
#include #define MAX 200 using namespace std; #define M_PI 3.1415926535897932384 int log2(int N) /*function to calculate the log2(.) of int numbers*/ { int k = N, i = 0; while(k) { k >>= 1; i++; } return i - 1; } int check(int n) //checking if the number of element is a power of 2 { return n > 0 && (n & (n - 1)) == 0; } int reverse(int N, int n) //calculating revers number { int j, p = 0; for(j = 1; j <= log2(N); j++) { if(n & (1 << (log2(N) - j))) p |= 1 << (j - 1); } return p; } void ordina(complex* f1, int N) //using the reverse order in the array { complex f2[MAX]; for(int i = 0; i < N; i++) f2[i] = f1[reverse(N, i)]; for(int j = 0; j < N; j++) f1[j] = f2[j]; } void transform(complex* f, int N) // { ordina(f, N); //first: reverse order complex *W; W = (complex *)malloc(N / 2 * sizeof(complex)); W[1] = polar(1., -2. * M_PI / N); W[0] = 1; for(int i = 2; i < N / 2; i++) W[i] = pow(W[1], i); int n = 1; int a = N / 2; for(int j = 0; j < log2(N); j++) { for(int i = 0; i < N; i++) { if(!(i & n)) { complex temp = f[i]; complex Temp = W[(i * a) % (n * a)] * f[i + n]; f[i] = temp + Temp; f[i + n] = temp - Temp; } } n *= 2; a = a / 2; } } void FFT(complex* f, int N, double d) { transform(f, N); for(int i = 0; i < N; i++) f[i] *= d; //multiplying by step } int main() { int n; do { cout << \"specify array dimension (MUST be power of 2)\" << endl; cin >> n; } while(!check(n)); double d; cout << \"specify sampling step\" << endl; //just write 1 in order to have the same results of matlab fft(.) cin >> d; complex vec[MAX]; cout << \"specify the array\" << endl; for(int i = 0; i < n; i++) { cout << \"specify element number: \" << i << endl; cin >> vec[i]; } FFT(vec, n, d); cout << \"...printing the FFT of the array specified\" << endl; for(int j = 0; j < n; j++) cout << vec[j] << endl; return 0; }