| /* |
| * 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; |
| } |