Tom Tsou | 3522229 | 2016-06-22 16:16:30 -0700 | [diff] [blame] | 1 | /* |
Pau Espin Pedrol | bdb970e | 2019-07-22 12:03:39 +0200 | [diff] [blame] | 2 | * Fast Fourier transform |
Tom Tsou | 3522229 | 2016-06-22 16:16:30 -0700 | [diff] [blame] | 3 | * |
| 4 | * Copyright (C) 2012 Tom Tsou <tom@tsou.cc> |
Pau Espin Pedrol | bdb970e | 2019-07-22 12:03:39 +0200 | [diff] [blame] | 5 | * |
Pau Espin Pedrol | 21d03d3 | 2019-07-22 12:05:52 +0200 | [diff] [blame^] | 6 | * SPDX-License-Identifier: LGPL-2.1+ |
| 7 | * |
Tom Tsou | 3522229 | 2016-06-22 16:16:30 -0700 | [diff] [blame] | 8 | * This program is free software: you can redistribute it and/or modify |
| 9 | * it under the terms of the GNU Affero General Public License as published by |
| 10 | * the Free Software Foundation, either version 3 of the License, or |
| 11 | * (at your option) any later version. |
Pau Espin Pedrol | bdb970e | 2019-07-22 12:03:39 +0200 | [diff] [blame] | 12 | * |
Tom Tsou | 3522229 | 2016-06-22 16:16:30 -0700 | [diff] [blame] | 13 | * This program is distributed in the hope that it will be useful, |
| 14 | * but WITHOUT ANY WARRANTY; without even the implied warranty of |
| 15 | * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the |
| 16 | * GNU Affero General Public License for more details. |
Pau Espin Pedrol | bdb970e | 2019-07-22 12:03:39 +0200 | [diff] [blame] | 17 | * |
Tom Tsou | 3522229 | 2016-06-22 16:16:30 -0700 | [diff] [blame] | 18 | * You should have received a copy of the GNU Affero General Public License |
| 19 | * along with this program; if not, see <http://www.gnu.org/licenses/>. |
| 20 | * See the COPYING file in the main directory for details. |
| 21 | */ |
| 22 | |
| 23 | #include <stdlib.h> |
| 24 | #include <string.h> |
| 25 | #include <assert.h> |
| 26 | #include <fftw3.h> |
| 27 | |
| 28 | #include "fft.h" |
| 29 | |
| 30 | struct fft_hdl { |
| 31 | float *fft_in; |
| 32 | float *fft_out; |
| 33 | int len; |
| 34 | fftwf_plan fft_plan; |
| 35 | }; |
| 36 | |
Pau Espin Pedrol | bdb970e | 2019-07-22 12:03:39 +0200 | [diff] [blame] | 37 | /*! \brief Initialize FFT backend |
Tom Tsou | 3522229 | 2016-06-22 16:16:30 -0700 | [diff] [blame] | 38 | * \param[in] reverse FFT direction |
Pau Espin Pedrol | bdb970e | 2019-07-22 12:03:39 +0200 | [diff] [blame] | 39 | * \param[in] m FFT length |
Tom Tsou | 3522229 | 2016-06-22 16:16:30 -0700 | [diff] [blame] | 40 | * \param[in] istride input stride count |
| 41 | * \param[in] ostride output stride count |
| 42 | * \param[in] in input buffer (FFTW aligned) |
| 43 | * \param[in] out output buffer (FFTW aligned) |
| 44 | * \param[in] ooffset initial offset into output buffer |
| 45 | * |
| 46 | * If the reverse is non-NULL, then an inverse FFT will be used. This is a |
| 47 | * wrapper for advanced non-contiguous FFTW usage. See FFTW documentation for |
| 48 | * further details. |
| 49 | * |
| 50 | * http://www.fftw.org/doc/Advanced-Complex-DFTs.html |
| 51 | * |
| 52 | * It is currently unknown how the offset of the output buffer affects FFTW |
| 53 | * memory alignment. |
| 54 | */ |
| 55 | struct fft_hdl *init_fft(int reverse, int m, int istride, int ostride, |
| 56 | float *in, float *out, int ooffset) |
| 57 | { |
| 58 | int rank = 1; |
| 59 | int n[] = { m }; |
| 60 | int howmany = istride; |
| 61 | int idist = 1; |
| 62 | int odist = 1; |
| 63 | int *inembed = n; |
| 64 | int *onembed = n; |
| 65 | fftwf_complex *obuffer, *ibuffer; |
| 66 | |
| 67 | struct fft_hdl *hdl = (struct fft_hdl *) malloc(sizeof(struct fft_hdl)); |
| 68 | if (!hdl) |
| 69 | return NULL; |
| 70 | |
| 71 | int direction = FFTW_FORWARD; |
| 72 | if (reverse) |
| 73 | direction = FFTW_BACKWARD; |
| 74 | |
| 75 | ibuffer = (fftwf_complex *) in; |
| 76 | obuffer = (fftwf_complex *) out + ooffset; |
| 77 | |
| 78 | hdl->fft_in = in; |
| 79 | hdl->fft_out = out; |
| 80 | hdl->fft_plan = fftwf_plan_many_dft(rank, n, howmany, |
| 81 | ibuffer, inembed, istride, idist, |
| 82 | obuffer, onembed, ostride, odist, |
| 83 | direction, FFTW_MEASURE); |
| 84 | return hdl; |
| 85 | } |
| 86 | |
| 87 | void *fft_malloc(size_t size) |
| 88 | { |
| 89 | return fftwf_malloc(size); |
| 90 | } |
| 91 | |
| 92 | void fft_free(void *ptr) |
| 93 | { |
| 94 | free(ptr); |
| 95 | } |
| 96 | |
Pau Espin Pedrol | bdb970e | 2019-07-22 12:03:39 +0200 | [diff] [blame] | 97 | /*! \brief Free FFT backend resources |
Tom Tsou | 3522229 | 2016-06-22 16:16:30 -0700 | [diff] [blame] | 98 | */ |
| 99 | void free_fft(struct fft_hdl *hdl) |
| 100 | { |
| 101 | fftwf_destroy_plan(hdl->fft_plan); |
| 102 | free(hdl); |
| 103 | } |
| 104 | |
| 105 | /*! \brief Run multiple DFT operations with the initialized plan |
| 106 | * \param[in] hdl handle to an intitialized fft struct |
| 107 | * |
| 108 | * Input and output buffers are configured with init_fft(). |
| 109 | */ |
| 110 | int cxvec_fft(struct fft_hdl *hdl) |
| 111 | { |
| 112 | fftwf_execute(hdl->fft_plan); |
| 113 | return 0; |
| 114 | } |