113 lines
2.9 KiB
C
113 lines
2.9 KiB
C
/*
|
|
* Fast Fourier transform
|
|
*
|
|
* Copyright (C) 2012 Tom Tsou <tom@tsou.cc>
|
|
*
|
|
* This program is free software: you can redistribute it and/or modify
|
|
* it under the terms of the GNU Affero General Public License as published by
|
|
* the Free Software Foundation, either version 3 of the License, or
|
|
* (at your option) any later version.
|
|
*
|
|
* This program is distributed in the hope that it will be useful,
|
|
* but WITHOUT ANY WARRANTY; without even the implied warranty of
|
|
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
|
|
* GNU Affero General Public License for more details.
|
|
*
|
|
* You should have received a copy of the GNU Affero General Public License
|
|
* along with this program; if not, see <http://www.gnu.org/licenses/>.
|
|
* See the COPYING file in the main directory for details.
|
|
*/
|
|
|
|
#include <stdlib.h>
|
|
#include <string.h>
|
|
#include <assert.h>
|
|
#include <fftw3.h>
|
|
|
|
#include "fft.h"
|
|
|
|
struct fft_hdl {
|
|
float *fft_in;
|
|
float *fft_out;
|
|
int len;
|
|
fftwf_plan fft_plan;
|
|
};
|
|
|
|
/*! \brief Initialize FFT backend
|
|
* \param[in] reverse FFT direction
|
|
* \param[in] m FFT length
|
|
* \param[in] istride input stride count
|
|
* \param[in] ostride output stride count
|
|
* \param[in] in input buffer (FFTW aligned)
|
|
* \param[in] out output buffer (FFTW aligned)
|
|
* \param[in] ooffset initial offset into output buffer
|
|
*
|
|
* If the reverse is non-NULL, then an inverse FFT will be used. This is a
|
|
* wrapper for advanced non-contiguous FFTW usage. See FFTW documentation for
|
|
* further details.
|
|
*
|
|
* http://www.fftw.org/doc/Advanced-Complex-DFTs.html
|
|
*
|
|
* It is currently unknown how the offset of the output buffer affects FFTW
|
|
* memory alignment.
|
|
*/
|
|
struct fft_hdl *init_fft(int reverse, int m, int istride, int ostride,
|
|
float *in, float *out, int ooffset)
|
|
{
|
|
int rank = 1;
|
|
int n[] = { m };
|
|
int howmany = istride;
|
|
int idist = 1;
|
|
int odist = 1;
|
|
int *inembed = n;
|
|
int *onembed = n;
|
|
fftwf_complex *obuffer, *ibuffer;
|
|
|
|
struct fft_hdl *hdl = (struct fft_hdl *) malloc(sizeof(struct fft_hdl));
|
|
if (!hdl)
|
|
return NULL;
|
|
|
|
int direction = FFTW_FORWARD;
|
|
if (reverse)
|
|
direction = FFTW_BACKWARD;
|
|
|
|
ibuffer = (fftwf_complex *) in;
|
|
obuffer = (fftwf_complex *) out + ooffset;
|
|
|
|
hdl->fft_in = in;
|
|
hdl->fft_out = out;
|
|
hdl->fft_plan = fftwf_plan_many_dft(rank, n, howmany,
|
|
ibuffer, inembed, istride, idist,
|
|
obuffer, onembed, ostride, odist,
|
|
direction, FFTW_MEASURE);
|
|
return hdl;
|
|
}
|
|
|
|
void *fft_malloc(size_t size)
|
|
{
|
|
return fftwf_malloc(size);
|
|
}
|
|
|
|
void fft_free(void *ptr)
|
|
{
|
|
free(ptr);
|
|
}
|
|
|
|
/*! \brief Free FFT backend resources
|
|
*/
|
|
void free_fft(struct fft_hdl *hdl)
|
|
{
|
|
fftwf_destroy_plan(hdl->fft_plan);
|
|
free(hdl);
|
|
}
|
|
|
|
/*! \brief Run multiple DFT operations with the initialized plan
|
|
* \param[in] hdl handle to an intitialized fft struct
|
|
*
|
|
* Input and output buffers are configured with init_fft().
|
|
*/
|
|
int cxvec_fft(struct fft_hdl *hdl)
|
|
{
|
|
fftwf_execute(hdl->fft_plan);
|
|
return 0;
|
|
}
|