From 008876dbabffaf737f081e03530656127eafe0c0 Mon Sep 17 00:00:00 2001 From: "J.A. de Jong" Date: Sat, 17 Feb 2018 17:01:53 +0100 Subject: [PATCH] Added inverse FFT code. --- beamforming/c/fft.c | 71 ++++++++++++++++++++++++++++++++++++++-- beamforming/c/fft.h | 15 +++++++-- beamforming/wrappers.pyx | 31 ++++++++++++++++++ test/fft_test.py | 13 +++++--- 4 files changed, 121 insertions(+), 9 deletions(-) diff --git a/beamforming/c/fft.c b/beamforming/c/fft.c index 870ae96..99478ee 100644 --- a/beamforming/c/fft.c +++ b/beamforming/c/fft.c @@ -50,6 +50,69 @@ void Fft_free(Fft* fft) { feTRACE(15); } us Fft_nfft(const Fft* fft) {return fft->nfft;} +void Fft_ifft_single(const Fft* fft,const vc* freqdata,vd* result) { + fsTRACE(15); + dbgassert(fft && freqdata && result,NULLPTRDEREF); + const us nfft = fft->nfft; + dbgassert(result->size == nfft, + "Invalid size for time data rows." + " Should be equal to nfft"); + + dbgassert(freqdata->size == (nfft/2+1),"Invalid number of rows in" + " result array"); + + + /* Obtain fft_result */ + vd fft_result = fft->fft_result; + + /* Copy freqdata, to fft_result. */ + d* fft_result_ptr = getvdval(&fft_result,0); + *fft_result_ptr = c_real(*getvcval(freqdata,0)); + + d_copy(&fft_result_ptr[1], + (d*) getvcval(freqdata,1), + nfft-1); + + /* Perform backward transform */ + npy_rfftb(nfft, + fft_result_ptr, + getvdval(&fft->fft_work,0)); + + /* Scale by dividing by nfft. Checked with numpy implementation + * that this indeed needs to be done. */ + d_scale(fft_result_ptr,1/((d) nfft),nfft); + + vd_copy(result, + &fft_result); + + feTRACE(15); +} +void Fft_ifft(const Fft* fft,const cmat* freqdata,dmat* timedata) { + fsTRACE(15); + + dbgassert(fft && timedata && freqdata,NULLPTRDEREF); + + const us nchannels = timedata->n_cols; + dbgassert(timedata->n_cols == freqdata->n_cols, + "Number of columns in timedata and result" + " should be equal."); + + for(us col=0;colfft_work); diff --git a/beamforming/c/fft.h b/beamforming/c/fft.h index 33bcf1e..e4073e7 100644 --- a/beamforming/c/fft.h +++ b/beamforming/c/fft.h @@ -57,19 +57,28 @@ void Fft_fft_single(const Fft* fft,const vd* timedata,vc* result); */ void Fft_fft(const Fft* fft,const dmat* timedata,cmat* result); +/** + * Perform inverse fft on a single channel. + * + * @param[in] fft Fft handle. + * @param[in] freqdata Frequency domain input data, to be iFft'th. + * @param[out] result: iFft't data, should have size (nfft). + */ +void Fft_ifft_single(const Fft* fft,const vc* freqdata,vd* result); + /** * Perform inverse FFT * * @param[in] fft Fft handle * @param[in] freqdata Frequency domain data - * @param[out] timedata + * @param[out] timedata Time domain result */ void Fft_ifft(const Fft* fft,const cmat* freqdata,dmat* timedata); /** - * + * Free up resources of Fft handle. * - * @param fft + * @param fft Fft handle. */ void Fft_free(Fft* fft); diff --git a/beamforming/wrappers.pyx b/beamforming/wrappers.pyx index 76a7ca8..63e90f6 100644 --- a/beamforming/wrappers.pyx +++ b/beamforming/wrappers.pyx @@ -20,6 +20,7 @@ cdef extern from "fft.h": c_Fft* Fft_alloc(us nfft) void Fft_free(c_Fft*) void Fft_fft(c_Fft*,dmat * timedate,cmat * res) nogil + void Fft_ifft(c_Fft*,cmat * freqdata,dmat* timedata) nogil us Fft_nfft(c_Fft*) @@ -64,6 +65,36 @@ cdef class Fft: return result + def ifft(self,c[::1,:] freqdata): + + cdef us nfft = Fft_nfft(self._fft) + cdef us nchannels = freqdata.shape[1] + assert freqdata.shape[0] == nfft//2+1 + + + # result[:,:] = np.nan+1j*np.nan + + cdef cmat f = cmat_foreign(freqdata.shape[0], + freqdata.shape[1], + &freqdata[0,0]) + + timedata = np.empty((nfft,nchannels), + dtype=NUMPY_FLOAT_TYPE, + order='F') + + cdef d[::1,:] timedata_view = timedata + cdef dmat t = dmat_foreign(timedata.shape[0], + timedata.shape[1], + &timedata_view[0,0]) + + Fft_ifft(self._fft,&f,&t) + + dmat_free(&t) + cmat_free(&f) + + return timedata + + cdef extern from "window.h": ctypedef enum WindowType: diff --git a/test/fft_test.py b/test/fft_test.py index 42a575d..0e7b041 100644 --- a/test/fft_test.py +++ b/test/fft_test.py @@ -8,14 +8,16 @@ Created on Mon Jan 15 19:45:33 2018 import numpy as np from beamforming import Fft -nfft=6 +nfft=9 print('nfft:',nfft) print(nfft) -nchannels = 1 +nchannels = 4 t = np.linspace(0,1,nfft+1)[:-1] # print(t) -x1 = 1+np.sin(2*np.pi*t)+3.2*np.cos(2*np.pi*t)+np.sin(7*np.pi*t) +#x1 = 1+np.sin(2*np.pi*t)+3.2*np.cos(2*np.pi*t)+np.sin(7*np.pi*t) +#x1 = np.sin(2*np.pi*t) +x1 = 1+0*t x = np.vstack([x1.T]*nchannels).T # Using transpose to get the strides right x = np.random.randn(nchannels,nfft).T @@ -28,6 +30,9 @@ print(X) fft = Fft(nfft) Y = fft.fft(x) -print('Fftpack fft') +print('Beamforming fft') print(Y) + +x2 = fft.ifft(Y) +print('normdiff:',np.linalg.norm(x2-x)) print('end python script')