Matlab code only! Write MATLAB or C code to perform a 28-point (N-28) Cooley-Tuk
ID: 2248976 • Letter: M
Question
Matlab code only!
Write MATLAB or C code to perform a 28-point (N-28) Cooley-Tukey FFT. Use a 4x7 array for the transform. Transforms on the individual rows and columns should be performed by direct implementation of the DFT as nested FOR loops. Compare the number of operations required for the FFT versus the DFT in this case. In your answer you should describe the steps required to perform the transform and explain the indexing of the data. Make sure you check that the FFT gives the same result as the direct application of the DFT (6 marks) Re-implement Question 5 using the Good-Thomas (Prime Factor) algorithm to perform the FFT. Comment on the performance of this implementation as compared to the Cooley- Tukey Method. (6 Marks)Explanation / Answer
#include <complex>
#include <iostream>
#include <valarray>
const double PI = 3.141592653589793238460;
typedef std::complex<double> Complex;
typedef std::valarray<Complex> CArray;
// Cooley–Tukey FFT
void fft(CArray& x)
{
const size_t N = x.size();
if (N <= 1) return;
// divide
CArray even = x[std::slice(0, N/2, 2)];
CArray odd = x[std::slice(1, N/2, 2)];
// conquer
fft(even);
fft(odd);
// combine
for (size_t k = 0; k < N/2; ++k)
{
Complex t = std::polar(1.0, -2 * PI * k / N) * odd[k];
x[k ] = even[k] + t;
x[k+N/2] = even[k] - t;
}
}
// Cooley-Tukey FFT (in-place, breadth-first, decimation-in-frequency)
void fft(CArray &x)
{
// DFT
unsigned int N = x.size(), k = N, n;
double thetaT = 3.14159265358979323846264338328L / N;
Complex phiT = Complex(cos(thetaT), -sin(thetaT)), T;
while (k > 1)
{
n = k;
k >>= 1;
phiT = phiT * phiT;
T = 1.0L;
for (unsigned int l = 0; l < k; l++)
{
for (unsigned int a = l; a < N; a += n)
{
unsigned int b = a + k;
Complex t = x[a] - x[b];
x[a] += x[b];
x[b] = t * T;
}
T *= phiT;
}
}
// Decimate
unsigned int m = (unsigned int)log2(N);
for (unsigned int a = 0; a < N; a++)
{
unsigned int b = a;
// Reverse bits
b = (((b & 0xaaaaaaaa) >> 1) | ((b & 0x55555555) << 1));
b = (((b & 0xcccccccc) >> 2) | ((b & 0x33333333) << 2));
b = (((b & 0xf0f0f0f0) >> 4) | ((b & 0x0f0f0f0f) << 4));
b = (((b & 0xff00ff00) >> 8) | ((b & 0x00ff00ff) << 8));
b = ((b >> 16) | (b << 16)) >> (32 - m);
if (b > a)
{
Complex t = x[a];
x[a] = x[b];
x[b] = t;
}
}
//// Normalize (This section make it not working correctly)
//Complex f = 1.0 / sqrt(N);
//for (unsigned int i = 0; i < N; i++)
// x[i] *= f;
}
// inverse fft (in-place)
void ifft(CArray& x)
{
// conjugate the complex numbers
x = x.apply(std::conj);
// forward fft
fft( x );
// conjugate the complex numbers again
x = x.apply(std::conj);
// scale the numbers
x /= x.size();
}
int main()
{
const Complex test[] = { 1.0, 1.0, 1.0, 1.0, 0.0, 0.0, 0.0, 0.0 };
CArray data(test, 28);
// forward fft
fft(data);
std::cout << "fft" << std::endl;
for (int i = 0; i < 28; ++i)
{
std::cout << data[i] << std::endl;
}
// inverse fft
ifft(data);
std::cout << std::endl << "ifft" << std::endl;
for (int i = 0; i < 28; ++i)
{
std::cout << data[i] << std::endl;
}
return 0;
}