I am trying to compute the FFT and then the IFFT just to try out if I can get the same signal back but I am not really sure how to accomplish it. This is how I do the FFT:
plan = fftw_plan_r2r_1d(blockSize, datas, out, FFTW_R2HC, FFTW_ESTIMATE);
fftw_execute(plan);
FFT (Fast Fourier Transform) is able to convert a signal from the time domain to the frequency domain. IFFT (Inverse FFT) converts a signal from the frequency domain to the time domain. The FFT of a non-periodic signal will cause the resulting frequency spectrum to suffer from leakage.
We believe that FFTW, which is free software, should become the FFT library of choice for most applications.
Code generation with MATLAB Coder™ supports fftw only for MEX output. For standalone C/C++ code, to select a planning method for FFT library calls, implement a getPlanMethod method in an FFT library callback class.
FFTW is a C subroutine library for computing the discrete Fourier transform (DFT) in one or more dimensions, of arbitrary input size, and of both real and complex data (as well as of even/odd data, i.e. the discrete cosine/sine transforms or DCT/DST).
Here is an example. It does two things. First, it prepares an input array in[N]
as a cosine wave, whose frequency is 3 and magnitude is 1.0, and Fourier transforms it. So, in the output, you should see a peak at out[3]
and and another at out[N-3]
. Since the magnitude of the cosine wave is 1.0, you get N/2 at out[3]
and out[N-3]
.
Second, it inverse Fourier transforms the array out[N]
back to in2[N]
. And after a proper normalization, you can see that in2[N]
is identical to in[N]
.
#include <stdlib.h>
#include <math.h>
#include <fftw3.h>
#define N 16
int main(void) {
fftw_complex in[N], out[N], in2[N]; /* double [2] */
fftw_plan p, q;
int i;
/* prepare a cosine wave */
for (i = 0; i < N; i++) {
in[i][0] = cos(3 * 2*M_PI*i/N);
in[i][1] = 0;
}
/* forward Fourier transform, save the result in 'out' */
p = fftw_plan_dft_1d(N, in, out, FFTW_FORWARD, FFTW_ESTIMATE);
fftw_execute(p);
for (i = 0; i < N; i++)
printf("freq: %3d %+9.5f %+9.5f I\n", i, out[i][0], out[i][1]);
fftw_destroy_plan(p);
/* backward Fourier transform, save the result in 'in2' */
printf("\nInverse transform:\n");
q = fftw_plan_dft_1d(N, out, in2, FFTW_BACKWARD, FFTW_ESTIMATE);
fftw_execute(q);
/* normalize */
for (i = 0; i < N; i++) {
in2[i][0] *= 1./N;
in2[i][1] *= 1./N;
}
for (i = 0; i < N; i++)
printf("recover: %3d %+9.5f %+9.5f I vs. %+9.5f %+9.5f I\n",
i, in[i][0], in[i][1], in2[i][0], in2[i][1]);
fftw_destroy_plan(q);
fftw_cleanup();
return 0;
}
This is the output:
freq: 0 -0.00000 +0.00000 I
freq: 1 +0.00000 +0.00000 I
freq: 2 -0.00000 +0.00000 I
freq: 3 +8.00000 -0.00000 I
freq: 4 +0.00000 +0.00000 I
freq: 5 -0.00000 +0.00000 I
freq: 6 +0.00000 -0.00000 I
freq: 7 -0.00000 +0.00000 I
freq: 8 +0.00000 +0.00000 I
freq: 9 -0.00000 -0.00000 I
freq: 10 +0.00000 +0.00000 I
freq: 11 -0.00000 -0.00000 I
freq: 12 +0.00000 -0.00000 I
freq: 13 +8.00000 +0.00000 I
freq: 14 -0.00000 -0.00000 I
freq: 15 +0.00000 -0.00000 I
Inverse transform:
recover: 0 +1.00000 +0.00000 I vs. +1.00000 +0.00000 I
recover: 1 +0.38268 +0.00000 I vs. +0.38268 +0.00000 I
recover: 2 -0.70711 +0.00000 I vs. -0.70711 +0.00000 I
recover: 3 -0.92388 +0.00000 I vs. -0.92388 +0.00000 I
recover: 4 -0.00000 +0.00000 I vs. -0.00000 +0.00000 I
recover: 5 +0.92388 +0.00000 I vs. +0.92388 +0.00000 I
recover: 6 +0.70711 +0.00000 I vs. +0.70711 +0.00000 I
recover: 7 -0.38268 +0.00000 I vs. -0.38268 +0.00000 I
recover: 8 -1.00000 +0.00000 I vs. -1.00000 +0.00000 I
recover: 9 -0.38268 +0.00000 I vs. -0.38268 +0.00000 I
recover: 10 +0.70711 +0.00000 I vs. +0.70711 +0.00000 I
recover: 11 +0.92388 +0.00000 I vs. +0.92388 +0.00000 I
recover: 12 +0.00000 +0.00000 I vs. +0.00000 +0.00000 I
recover: 13 -0.92388 +0.00000 I vs. -0.92388 +0.00000 I
recover: 14 -0.70711 +0.00000 I vs. -0.70711 +0.00000 I
recover: 15 +0.38268 +0.00000 I vs. +0.38268 +0.00000 I
Here's what I made. Mine is designed specifically to give the same output as Matlab. This only works on a column matrix (you can replace AMatrix
with a std::vector
).
AMatrix<complex<double>> FastFourierTransform::Forward1d(const AMatrix<double> &inMat)
{
AMatrix<complex<double>> outMat{ inMat.NumRows(), inMat.NumCols() };
size_t N = inMat.NumElements();
bool isOdd = N % 2 == 1;
size_t outSize = (isOdd) ? ceil(N / 2 + 1) : N / 2;
fftw_plan plan = fftw_plan_dft_r2c_1d(
inMat.NumRows(),
reinterpret_cast<double*>(&inMat.mutData()[0]), // mutData here is as same v.data() for std::vector
reinterpret_cast<fftw_complex*>(&outMat.mutData()[0]),
FFTW_ESTIMATE);
fftw_execute(plan);
// mirror
auto halfWayPoint = outMat.begin() + (outMat.NumRows() / 2) + 1;
auto startCopyDest = (isOdd) ? halfWayPoint : halfWayPoint - 1;
std::reverse_copy(outMat.begin() + 1, halfWayPoint, startCopyDest); // skip DC (+1)
std::for_each(halfWayPoint, outMat.end(), [](auto &c) { return c = std::conj(c);});
return std::move(outMat);
}
AMatrix<complex<double>> FastFourierTransform::Forward1d(const AMatrix<double> &inMat, size_t points)
{
// append zeros to matrix until it's the same size as points
AMatrix<double> sig = inMat;
sig.Resize(points, sig.NumCols());
for (size_t i = inMat.NumRows(); i < points; i++)
{
sig(i, 0) = 0;
}
return Forward1d(sig);
}
AMatrix<complex<double>> FastFourierTransform::Inverse1d(const AMatrix<complex<double>> &inMat)
{
AMatrix<complex<double>> outMat{ inMat.NumRows(), inMat.NumCols() };
size_t N = inMat.NumElements();
fftw_plan plan = fftw_plan_dft_1d(
inMat.NumRows(),
reinterpret_cast<fftw_complex*>(&inMat.mutData()[0]),
reinterpret_cast<fftw_complex*>(&outMat.mutData()[0]),
FFTW_BACKWARD,
FFTW_ESTIMATE);
fftw_execute(plan);
// Matlab normalizes
auto normalize = [=](auto &c) { return c *= 1. / N; };
std::for_each(outMat.begin(), outMat.end(), normalize);
fftw_destroy_plan(plan);
return std::move(outMat);
}
// From Matlab documentation: "ifft tests X to see whether vectors in X along the active dimension
// are conjugate symmetric. If so, the computation is faster and the output is real.
// An N-element vector x is conjugate symmetric if x(i) = conj(x(mod(N-i+1,N)+1)) for each element of x."
// http://uk.mathworks.com/help/matlab/ref/ifft.html
AMatrix<double> FastFourierTransform::Inverse1dConjSym(const AMatrix<complex<double>> &inMat)
{
AMatrix<double> outMat{ inMat.NumRows(), inMat.NumCols() };
size_t N = inMat.NumElements();
fftw_plan plan = fftw_plan_dft_c2r_1d(
inMat.NumRows(),
reinterpret_cast<fftw_complex*>(&inMat.mutData()[0]),
reinterpret_cast<double*>(&outMat.mutData()[0]),
FFTW_BACKWARD | FFTW_ESTIMATE);
fftw_execute(plan);
auto normalize = [=](auto &c) { return c *= 1. / N; };
std::for_each(outMat.begin(), outMat.end(), normalize);
fftw_destroy_plan(plan);
return std::move(outMat);
}
If you love us? You can donate to us via Paypal or buy me a coffee so we can maintain and grow! Thank you!
Donate Us With