diff --git a/.gitmodules b/.gitmodules index 61100c0..89efc2b 100644 --- a/.gitmodules +++ b/.gitmodules @@ -16,3 +16,6 @@ [submodule "third_party/lockfree"] path = third_party/lockfree url = https://github.com/boostorg/lockfree +[submodule "third_party/carma"] + path = third_party/carma + url = https://github.com/RUrlus/carma diff --git a/CMakeLists.txt b/CMakeLists.txt index c447cdf..74ef5a5 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -26,7 +26,7 @@ option(LASP_HAS_RTAUDIO "Compile with RtAudio Daq backend" ON) option(LASP_HAS_ULDAQ "Compile with UlDaq backend" ON) set(LASP_FFT_BACKEND "FFTW" CACHE STRING "FFT Library backend") -set_property(CACHE LASP_FFT_BACKEND PROPERTY STRINGS "FFTW" "FFTPack" "None") +set_property(CACHE LASP_FFT_BACKEND PROPERTY STRINGS "FFTW" "Armadillo") set(PY_FULL_VERSION ${PROJECT_VERSION}${PY_VERSION_SUFFIX}) @@ -52,10 +52,16 @@ include(OSSpecific) set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -Wall -Wextra -Wno-type-limits") -set(CMAKE_C_FLAGS_RELEASE "-O3 -mfpmath=sse -march=x86-64 -mtune=native \ +set(CMAKE_C_FLAGS_RELEASE "-O3 -flto -mfpmath=sse -march=x86-64 -mtune=native \ -fdata-sections -ffunction-sections -fomit-frame-pointer -finline-functions") # ############################# End compilation flags +add_subdirectory(third_party/carma) -add_subdirectory(third_party) +if(LASP_FFT_BACKEND STREQUAL "FFTW") + find_library(FFTW_LIBRARY NAMES fftw3 fftw) + set(FFT_LIBRARIES "${FFTW_LIBRARY}") +elseif(LASP_FFT_BACKEND STREQUAL "Armadillo") +endif() add_subdirectory(src/lasp) + diff --git a/src/lasp/__init__.py b/src/lasp/__init__.py index 28151ff..a101f99 100644 --- a/src/lasp/__init__.py +++ b/src/lasp/__init__.py @@ -16,7 +16,7 @@ from .lasp_measurement import * # Measurement, scaleBlockSens from .lasp_octavefilter import * # from .lasp_slm import * # SLM, Dummy from .lasp_record import * # RecordStatus, Recording -from .lasp_daqconfigs import DaqConfigurations +from .lasp_daqconfigs import * # from .lasp_siggen import * # SignalType, NoiseType, SiggenMessage, SiggenData, Siggen # from .lasp_weighcal import * # WeighCal # from .tools import * # SmoothingType, smoothSpectralData, SmoothingWidth diff --git a/src/lasp/c/lasp_fft.c b/src/lasp/c/lasp_fft.c deleted file mode 100644 index 310897d..0000000 --- a/src/lasp/c/lasp_fft.c +++ /dev/null @@ -1,263 +0,0 @@ -// lasp_fft.c -// -// Author: J.A. de Jong - ASCEE -// -// Description: -// -////////////////////////////////////////////////////////////////////// -#define TRACERPLUS (-5) -#include "lasp_types.h" -#include "lasp_config.h" -#include "lasp_tracer.h" -#include "lasp_fft.h" - -#if LASP_FFT_BACKEND == FFTPack -#include "fftpack.h" -typedef struct Fft_s { - us nfft; - vd fft_work; // Storage memory for fftpack -} Fft_s; -#elif LASP_FFT_BACKEND == FFTW -#include - -typedef struct Fft_s { - us nfft; - fftw_plan forward_plan; - fftw_plan reverse_plan; - c* complex_storage; - d* real_storage; -} Fft_s; -#else -#error "Cannot compile lasp_ffc.c, no FFT backend specified. Should either be FFTPack, or FFTW" -#endif - -void load_fft_wisdom(const char* wisdom) { -#if LASP_FFT_BACKEND == FFTPack -#elif LASP_FFT_BACKEND == FFTW - if(wisdom) { - int rv= fftw_import_wisdom_from_string(wisdom); - if(rv != 1) { - fprintf(stderr, "Error loading FFTW wisdom"); - } - } -#endif -} - -char* store_fft_wisdom() { -#if LASP_FFT_BACKEND == FFTPack - return NULL; -#elif LASP_FFT_BACKEND == FFTW - return fftw_export_wisdom_to_string(); -#endif -} - -Fft* Fft_create(const us nfft) { - fsTRACE(15); - if(nfft == 0) { - WARN("nfft should be > 0"); - return NULL; - } - - Fft* fft = a_malloc(sizeof(Fft)); - - fft->nfft = nfft; - -#if LASP_FFT_BACKEND == FFTPack - /* Initialize foreign fft lib */ - fft->fft_work = vd_alloc(2*nfft+15); - npy_rffti(nfft,getvdval(&fft->fft_work,0)); - check_overflow_vx(fft->fft_work); -#elif LASP_FFT_BACKEND == FFTW - fft->complex_storage = fftw_malloc(sizeof(c) * (nfft/2 + 1)); - fft->real_storage = fftw_malloc(sizeof(d) * nfft); - - fft->forward_plan = fftw_plan_dft_r2c_1d(nfft, - fft->real_storage, - fft->complex_storage, - FFTW_MEASURE); - fft->reverse_plan = fftw_plan_dft_c2r_1d(nfft, - fft->complex_storage, - fft->real_storage, - FFTW_MEASURE); - -#endif - - /* print_vd(&fft->fft_work); */ - - feTRACE(15); - return fft; -} -void Fft_free(Fft* fft) { - fsTRACE(15); - dbgassert(fft,NULLPTRDEREF); -#if LASP_FFT_BACKEND == FFTPack - vd_free(&fft->fft_work); -#elif LASP_FFT_BACKEND == FFTW - fftw_free(fft->complex_storage); - fftw_free(fft->real_storage); - fftw_destroy_plan(fft->forward_plan); - fftw_destroy_plan(fft->reverse_plan); -#endif - a_free(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->n_rows == nfft, - "Invalid size for time data rows." - " Should be equal to nfft"); - - dbgassert(freqdata->n_rows == (nfft/2+1),"Invalid number of rows in" - " result array"); - - - d* result_ptr = getvdval(result,0); - -#if LASP_FFT_BACKEND == FFTPack - d* freqdata_ptr = (d*) getvcval(freqdata,0); - /* Copy freqdata, to fft_result. */ - d_copy(&result_ptr[1],&freqdata_ptr[2],nfft-1,1,1); - result_ptr[0] = freqdata_ptr[0]; - - /* Perform inplace backward transform */ - npy_rfftb(nfft, - result_ptr, - getvdval(&fft->fft_work,0)); - - -#elif LASP_FFT_BACKEND == FFTW - c* freqdata_ptr = (c*) getvcval(freqdata,0); - - c_copy(fft->complex_storage, freqdata_ptr,nfft/2+1); - - fftw_execute(fft->reverse_plan); - - d_copy(result_ptr, fft->real_storage, nfft, 1, 1); - -#endif - check_overflow_vx(*result); - - /* Scale by dividing by nfft. Checked with numpy implementation - * that this indeed needs to be done for FFTpack. */ - d_scale(result_ptr,1/((d) nfft),nfft); - 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;colnfft; - assert_vx(timedata); - assert_vx(result); - dbgassert(timedata->n_rows == nfft, - "Invalid size for time data rows." - " Should be equal to nfft"); - - dbgassert(result->n_rows == (nfft/2+1),"Invalid number of rows in" - " result array"); - - -#if LASP_FFT_BACKEND == FFTPack - d* result_ptr = (d*) getvcval(result,0); - - /* Fftpack stores the data a bit strange, the resulting array - * has the DC value at 0,the first cosine at 1, the first sine - * at 2 etc. 1 - * resulting matrix, as for the complex data, the imaginary - * part of the DC component equals zero. */ - - /* Copy timedata, as it will be overwritten in the fft pass. */ - d_copy(&result_ptr[1],getvdval(timedata,0),nfft,1,1); - - - /* Perform fft */ - npy_rfftf(nfft,&result_ptr[1], - getvdval(&fft->fft_work,0)); - - /* Set real part of DC component to first index of the rfft - * routine */ - result_ptr[0] = result_ptr[1]; - - result_ptr[1] = 0; /* Set imaginary part of DC component - * to zero */ - - /* For an even fft, the imaginary part of the Nyquist frequency - * bin equals zero.*/ - if(islikely(nfft%2 == 0)) { - result_ptr[nfft+1] = 0; - } - check_overflow_vx(fft->fft_work); -#elif LASP_FFT_BACKEND == FFTW - - d* timedata_ptr = getvdval(timedata,0); - c* result_ptr = getvcval(result,0); - d_copy(fft->real_storage,timedata_ptr, nfft, 1, 1); - - fftw_execute(fft->forward_plan); - - c_copy(result_ptr, fft->complex_storage, nfft/2+1); - -#endif - check_overflow_vx(*result); - feTRACE(15); - -} -void Fft_fft(const Fft* fft,const dmat* timedata,cmat* result) { - fsTRACE(15); - - dbgassert(fft && timedata && result,NULLPTRDEREF); - - const us nchannels = timedata->n_cols; - dbgassert(timedata->n_cols == result->n_cols, - "Number of columns in timedata and result" - " should be equal."); - - for(us col=0;col -#include #include /** diff --git a/src/lasp/dsp/CMakeLists.txt b/src/lasp/dsp/CMakeLists.txt index b5d679c..d951e78 100644 --- a/src/lasp/dsp/CMakeLists.txt +++ b/src/lasp/dsp/CMakeLists.txt @@ -5,6 +5,8 @@ set(lasp_dsp_files lasp_siggen.cpp lasp_siggen_impl.cpp lasp_window.cpp + lasp_fft.cpp + lasp_avpowerspectra.cpp ) diff --git a/src/lasp/dsp/lasp_avpowerspectra.cpp b/src/lasp/dsp/lasp_avpowerspectra.cpp new file mode 100644 index 0000000..d98d19d --- /dev/null +++ b/src/lasp/dsp/lasp_avpowerspectra.cpp @@ -0,0 +1,68 @@ +#define DEBUGTRACE_ENABLED +#include "lasp_avpowerspectra.h" +#include "debugtrace.hpp" +#include + +using std::runtime_error; + +PowerSpectra::PowerSpectra(const us nfft, const Window::WindowType w) + : PowerSpectra(Window::create(w, nfft)) {} + +PowerSpectra::PowerSpectra(const vd &window) + : nfft(window.size()), _fft(nfft), _window(window) { + + d win_pow = arma::sum(window % window); + + /* Scale fft such that power is easily computed */ + _scale_fac = sqrt(2 / win_pow) / nfft; + } +arma::Cube PowerSpectra::compute(const dmat &input) { + + dmat input_tmp = input; + + // Multiply each column of the inputs element-wise with the window. + input_tmp.each_col() %= _window; + + cmat rfft = _fft.fft(input_tmp) * _scale_fac; + + arma::cx_cube output(rfft.n_rows, input.n_cols, input.n_cols); + for (us i = 0; i < input.n_cols; i++) { + for (us j = 0; j < input.n_cols; j++) { + output.slice(j).col(i) = rfft.col(i) % arma::conj(rfft.col(j)); + } + } + return output; +} + +AvPowerSpectra::AvPowerSpectra(const us nfft, const Window::WindowType w, + const d overlap_percentage, + const int time_constant) + : _ps(nfft, w) { + + DEBUGTRACE_ENTER; + if (overlap_percentage >= 100 || overlap_percentage < 0) { + throw runtime_error("Overlap percentage should be >= 0 and < 100"); + } + + _overlap_jump = nfft - (nfft * overlap_percentage) / 100; + if (_overlap_jump == 0) { + throw runtime_error("Overlap is too high. Results in no jump. Please " + "choose a smaller overlap percentage or a higher nfft"); + } + if(time_constant < 0) { + _mode = Mode::Averaging; + } else if(time_constant == 0) { + _mode = Mode::Spectrogram; + } + else { + _mode = Mode::Leaking; + _alpha = exp(-static_cast(_overlap_jump)/time_constant); + } + + } + +arma::cx_cube* AvPowerSpectra::compute(const dmat& timedata) { + + + return nullptr; +} diff --git a/src/lasp/dsp/lasp_avpowerspectra.h b/src/lasp/dsp/lasp_avpowerspectra.h new file mode 100644 index 0000000..67acd6f --- /dev/null +++ b/src/lasp/dsp/lasp_avpowerspectra.h @@ -0,0 +1,155 @@ +#pragma once +#include +#include "lasp_fft.h" +#include "lasp_mathtypes.h" +#include "lasp_window.h" +#include + +/** + * @brief Computes single-sided cross-power spectra for a group of channels + */ +class PowerSpectra { +public: + /** + * @brief The FFT length + */ + us nfft; + +private: + /** + * @brief Instance to compute actual FFT's + */ + Fft _fft; + + vd _window; + c _scale_fac; + +public: + /** + * @brief Initalize power spectra computer + * + * @param nfft The fft length + * @param w The window type + */ + PowerSpectra(const us nfft, const Window::WindowType w); + + /** + * @brief Initalize power spectra computer + * + * @param window Uses the window length to compute fft length, and uses the + * window shape as the actual window. + */ + PowerSpectra(const vd &window); + + /** + * @brief Computes the spectra. Data is normalized by the (spectral) power of + * the window, to compensate for the effect of the window. + * + * @param input Input data, first axis is assumed the sample, second axis the + * channel. Input first dimension should be equal to nfft, otherwise a + * runtime error is thrown. + * + * @return Cross-power-spectra. Cubic array with size indices C_fij, where + * f is the frequency index, i is the row and j is the column. An element + * can be accessed by: C(f, i, j). Storage is such that frequency components + * are contiguous. + */ + arma::Cube compute(const dmat &input); +}; + +/** + * @brief Estimate cross-power spectra using Welch' method of spectral + * estimation. The exact amount of overlap in Welch' method is rounded up to a + * certain amount of samples. + */ +class AvPowerSpectra { + + enum class Mode { + Averaging = 0, // Averaging all time date + Leaking = 1, // Exponential weighting of an "instantaneous cps" + Spectrogram = 2 // Instantenous spectrum, no averaging + }; + + Mode _mode; + d _alpha; // Only valid in case of 'Leaking' + + PowerSpectra _ps; + /** + * @brief Current estimate of cross-spectral density + */ + arma::cx_cube _est; + + /** + * @brief Buffer of storage of time data. + */ + dmat _timeBuf; + /** + * @brief Current time index in ring buffer + */ + us _curTimeIdx = 0; + /** + * @brief The amount of samples to jump in the overlap + */ + us _overlap_jump; + +public: + /** + * @brief Initalize averaged power spectra computer. If a time constant is + * given > 0, it is used in a kind of exponential weighting. + * + * @param nfft The fft length + * @param w The window type. + * @param overlap_percentage A number 0 < overlap_percentage <= 100. It + * determines the amount of overlap used in Welch' method. A typical value is + * 50 %, i.e. 50. + * @param time_constant Value should either be < 0, indicating that the + * estimate is averages over all time data. + * For a value = 0 the instantaneous power spectrum is returned, which can be + * interpreted as the spectrogram. For a value > 0 a exponential forgetting is + * used, where the value is used as the time constant such that the decay + * follows approximately the trend exp(-n/time_constant), where n is the + * sample number in the power spectra. To choose 'fast' time weighting, set + * time_constant to the value of fs*0.125, where fs denotes the sampling + * frequency. + **/ + AvPowerSpectra(const us nfft = 2048, + const Window::WindowType w = Window::WindowType::Hann, + const d overlap_percentage = 50., + const int time_constant = -1); + + AvPowerSpectra(const AvPowerSpectra &) = delete; + AvPowerSpectra &operator=(const AvPowerSpectra &) = delete; + + void reset() { _est.reset(); } + + /** + * @brief Comnpute an update of the power spectra based on given time data. + * Note that the number of channels is determined from the first time this + * function is called. If a later call has an incompatible number of + * channels, a runtime error is thrown. + * + * @param timedata + * + * @return Optionally, a reference (NOT OWNING POINTER) to a new estimate of + * the power spectra. An update is only given if the amount of new time data + * is enough to compute a new estimate. This is the case when the amount of + * new time data is enough. + */ + arma::cx_cube *compute(const dmat &timedata); + + /** + * @brief Returns the latest estimate of cps (cross-power spectra. + * + * @return Pointer (reference, not owning!) to spectral estimate date, or + * nullptr, in case nothing could yet be computed. + */ + arma::cx_cube *get_est(); + + /** + * @brief The overlap is rounded to a certain amount of time samples. This + * function returns that value. + * + * @return The amount of samples in overlapping. + */ + us exactOverlapSamples() const { return _ps.nfft - _overlap_jump; } +}; diff --git a/src/lasp/dsp/lasp_biquadbank.cpp b/src/lasp/dsp/lasp_biquadbank.cpp new file mode 100644 index 0000000..e69de29 diff --git a/src/lasp/dsp/lasp_biquadbank.h b/src/lasp/dsp/lasp_biquadbank.h new file mode 100644 index 0000000..e69de29 diff --git a/src/lasp/dsp/lasp_fft.cpp b/src/lasp/dsp/lasp_fft.cpp new file mode 100644 index 0000000..46bdef9 --- /dev/null +++ b/src/lasp/dsp/lasp_fft.cpp @@ -0,0 +1,151 @@ +// lasp_fft.c +// +// Author: J.A. de Jong - ASCEE +// +// Description: +// FFt implementation +// +////////////////////////////////////////////////////////////////////// +#include +#define DEBUGTRACE_ENABLED +#include "lasp_fft.h" +#include "debugtrace.hpp" +#include "lasp_config.h" +using std::runtime_error; + +#if LASP_FFT_BACKEND == Armadillo +#include "fftpack.h" +class Fft_impl { +public: + us nfft; + Fft_impl(const us nfft) : nfft(nfft) { + throw runtime_error( + "This code does not output correct results, as it computes the full " + "FFT. It needs to be reworked to single-sided spectra."); + + DEBUGTRACE_ENTER; + } + + vd ifft(const vc &f) { return arma::ifft(f); } + vc fft(const vd &time) { return arma::fft(time); } +}; +#elif LASP_FFT_BACKEND == FFTW +#include + +/** + * @brief FFTW implementation + */ +class Fft_impl { +public: + us nfft; + fftw_plan forward_plan = nullptr; + fftw_plan reverse_plan = nullptr; + fftw_complex *frequencyDomain = nullptr; + d *timeDomain = nullptr; + + Fft_impl(const us nfft) : nfft(nfft) { + + timeDomain = (d *)fftw_malloc(sizeof(d) * nfft); + + frequencyDomain = + (fftw_complex *)fftw_malloc(sizeof(fftw_complex) * (nfft / 2 + 1)); + + forward_plan = + fftw_plan_dft_r2c_1d(nfft, timeDomain, frequencyDomain, FFTW_MEASURE); + + reverse_plan = + fftw_plan_dft_c2r_1d(nfft, frequencyDomain, timeDomain, FFTW_MEASURE); + }; + + ~Fft_impl() { + fftw_destroy_plan(forward_plan); + fftw_destroy_plan(reverse_plan); + fftw_free(frequencyDomain); + fftw_free(timeDomain); + } + + vc fft(const vd &time) { + + vc res(nfft / 2 + 1); + memcpy(timeDomain, time.memptr(), sizeof(d) * nfft); + + fftw_execute(forward_plan); + + memcpy(reinterpret_cast(res.memptr()), frequencyDomain, + sizeof(d) * (nfft / 2 + 1)); + + return res; + } + vd ifft(const vc &f) { + vd res(nfft); + memcpy(frequencyDomain, + reinterpret_cast(const_cast(f.memptr())), + (nfft / 2 + 1) * sizeof(c)); + + fftw_execute(reverse_plan); + + memcpy(res.memptr(), frequencyDomain, nfft * sizeof(d)); + + /* Scale by dividing by nfft. Checked with numpy implementation + * that this indeed needs to be done for FFTW. */ + res *= 1 / nfft; + + return res; + } +}; +#else +#error \ + "Cannot compile lasp_ffc.c, no FFT backend specified. Should either be FFTPack, or FFTW" +#endif + +Fft::Fft(const us nfft) { + DEBUGTRACE_ENTER; + if (nfft == 0) { + throw std::runtime_error("Invalid nfft: 0"); + } + if (nfft >= LASP_MAX_NFFT) { + throw std::runtime_error("Invalid nfft, should be smaller than: " + + std::to_string(LASP_MAX_NFFT)); + } + _impl = std::make_unique(nfft); +} +Fft::~Fft() { } + +us Fft::nfft() const { return _impl->nfft; } + +vd Fft::ifft(const vc &freqdata) { + DEBUGTRACE_ENTER; + if (freqdata.size() != (_impl->nfft / 2 + 1)) { + throw runtime_error("Invalid size for frequency domain data rows." + " Should be equal to floor(nfft/2)+1"); + } + return _impl->ifft(freqdata); +} +dmat Fft::ifft(const cmat &freqdata) { + dmat res(_impl->nfft, freqdata.n_cols); + for (us colno = 0; colno < freqdata.n_cols; colno++) { + res.col(colno) = _impl->ifft(freqdata.col(colno)); + } + return res; +} +vc Fft::fft(const vd &timedata) { return _impl->fft(timedata); } + +void load_fft_wisdom(const char *wisdom) { +#if LASP_FFT_BACKEND == Armadillo +#elif LASP_FFT_BACKEND == FFTW + if (wisdom) { + int rv = fftw_import_wisdom_from_string(wisdom); + if (rv != 1) { + fprintf(stderr, "Error loading FFTW wisdom"); + } + } +#endif +} + +char *store_fft_wisdom() { +#if LASP_FFT_BACKEND == Armadillo + return NULL; +#elif LASP_FFT_BACKEND == FFTW + return fftw_export_wisdom_to_string(); +#endif +} diff --git a/src/lasp/dsp/lasp_fft.h b/src/lasp/dsp/lasp_fft.h new file mode 100644 index 0000000..7d99489 --- /dev/null +++ b/src/lasp/dsp/lasp_fft.h @@ -0,0 +1,90 @@ +// lasp_fft.h +// +// Author: J.A. de Jong - ASCEE +// +// Description: FFT class +// Interface to the FFT library, multiple channel FFT's +////////////////////////////////////////////////////////////////////// +#pragma once +#include +#include "lasp_mathtypes.h" + +/** + * Load wisdom data from a NULL-terminated character string. Note that at this + * point in time this is only used by the FFTW fft backend. + * + * @param[in] wisdom: Wisdom string + */ +void load_fft_wisdom(const char* wisdom); + +/** + * Creates a NULL-terminated string containing FFTW wisdom, or state knowledge + * from other libraries. + * + * @returns A NULL-terminated string. *Should be deallocated after being used* + */ +char* store_fft_wisdom(void); + +class Fft_impl; + +/** + * Perform forward FFT's on real time data. Computes single-sided spectra. + */ +class Fft { + std::unique_ptr _impl; + + public: + /** + * @brief Initialize FFT + * + * @param nfft The length of nfft + */ + Fft(const us nfft); + ~Fft(); + + Fft(const Fft&) = delete; + Fft& operator=(const Fft&) = delete; + + /** + * @brief Return nfft + * + * @return nfft + */ + us nfft() const; + + /** + * Compute the fft for a single channel of data. + * + * @param[in] timedata Input time data, should have size nfft + * @returns result Result complex vector + * */ + vc fft(const vd& timedata); + + /** + * Compute the fft of the data matrix, first axis is assumed to be + * the time axis. + * + * @param[in] timedata Input time data, should have size nfft. First axis + * is time, second axis is channel + * @returns Result complex array + */ + cmat fft(const dmat& timedata); + + /** + * Perform inverse fft on a single channel. + * + * @param[in] freqdata Frequency domain input data, to be iFft'th. Should + * have size nfft/2+1 + * @returns timedata: iFft't data, size nfft. + */ + vd ifft(const vc& freqdata); + + /** + * Perform inverse FFT + * + * @param[in] freqdata Frequency domain data + * @param[out] timedata Time domain result + */ + dmat ifft(const cmat& freqdata); + +}; diff --git a/src/lasp/dsp/lasp_mathtypes.h b/src/lasp/dsp/lasp_mathtypes.h index 6ecce1d..f7078e2 100644 --- a/src/lasp/dsp/lasp_mathtypes.h +++ b/src/lasp/dsp/lasp_mathtypes.h @@ -40,6 +40,8 @@ #endif // LASP_DOUBLE_PRECISION using vd = arma::Col; +using vc = arma::Col; using dmat = arma::Mat; +using cmat = arma::Mat; const d number_pi = arma::datum::pi; diff --git a/src/lasp/lasp_cpp.cpp b/src/lasp/lasp_cpp.cpp index 03df0b1..fc3ff3c 100644 --- a/src/lasp/lasp_cpp.cpp +++ b/src/lasp/lasp_cpp.cpp @@ -1,5 +1,23 @@ #include "lasp_config.h" #include + +/*! \mainpage + * + * \section intro_sec Introduction + * + * Welcome to the LASP (Library for Acoustic Signal Processing) code + * documentation. The code comprises a part which is written in C++, a part + * that is written in Python, and a part that functions as glue, which is + * Pybind11 C++ glue code. An example of such a file is the current one. + * + * \section Installation + * + * For the installation manual, please refer to the README.md of the Git + * repository. It is recommended to install the software from source. + * + * + * */ + namespace py = pybind11; void init_dsp(py::module &m); @@ -18,9 +36,12 @@ PYBIND11_MODULE(lasp_cpp, m) { init_streammgr(m); init_datahandler(m); + // We store the version number of the code via CMake, and create an + // attribute in the C++ code. m.attr("__version__") = std::to_string(LASP_VERSION_MAJOR) + "." + std::to_string(LASP_VERSION_MINOR); + m.attr("LASP_VERSION_MAJOR") = LASP_VERSION_MAJOR; m.attr("LASP_VERSION_MINOR") = LASP_VERSION_MINOR; } diff --git a/third_party/CMakeLists.txt b/third_party/CMakeLists.txt deleted file mode 100644 index 08bfc3d..0000000 --- a/third_party/CMakeLists.txt +++ /dev/null @@ -1,10 +0,0 @@ -# third_party/CMakeLists.txt - -if(LASP_FFT_BACKEND STREQUAL "FFTW") - find_library(FFTW_LIBRARY NAMES fftw3 fftw) - set(FFT_LIBRARIES "${FFTW_LIBRARY}") -elseif(LASP_FFT_BACKEND STREQUAL "FFTPack") - include_directories(fftpack) - set(FFT_LIBRARIES fftpack) - add_subdirectory(fftpack) -endif() diff --git a/third_party/carma b/third_party/carma new file mode 160000 index 0000000..5673bc7 --- /dev/null +++ b/third_party/carma @@ -0,0 +1 @@ +Subproject commit 5673bc713e9515a00fba0c7378375fc075bdc693 diff --git a/third_party/fftpack/CMakeLists.txt b/third_party/fftpack/CMakeLists.txt deleted file mode 100644 index b9b9903..0000000 --- a/third_party/fftpack/CMakeLists.txt +++ /dev/null @@ -1,9 +0,0 @@ -# We borrow Numpy's implementation for doing the Fast Fourier Transform. -# This FFT code appears to be faster than KISSFFT. -add_library(fftpack - fftpack.c -) - -# Ling fft to math -target_link_libraries(fftpack PRIVATE m) -target_include_directories(fftpack PUBLIC .) diff --git a/third_party/fftpack/fftpack.c b/third_party/fftpack/fftpack.c deleted file mode 100644 index f5fe21f..0000000 --- a/third_party/fftpack/fftpack.c +++ /dev/null @@ -1,1501 +0,0 @@ -/* - * fftpack.c : A set of FFT routines in C. - * Algorithmically based on Fortran-77 FFTPACK by Paul N. Swarztrauber (Version 4, 1985). -*/ -#define NPY_VISIBILITY_HIDDEN -/* #define NPY_NO_DEPRECATED_API NPY_API_VERSION */ - -#include -#include - - -#define DOUBLE -#ifdef DOUBLE -#define Treal double -#else -#define Treal float -#endif - - -#define ref(u,a) u[a] - -#define MAXFAC 13 /* maximum number of factors in factorization of n */ -#define NSPECIAL 4 /* number of factors for which we have special-case routines */ - -#ifdef __cplusplus -extern "C" { -#endif - - -/* ---------------------------------------------------------------------- - passf2, passf3, passf4, passf5, passf. Complex FFT passes fwd and bwd. ------------------------------------------------------------------------ */ - -static void passf2(int ido, int l1, const Treal cc[], Treal ch[], const Treal wa1[], int isign) - /* isign==+1 for backward transform */ - { - int i, k, ah, ac; - Treal ti2, tr2; - if (ido <= 2) { - for (k=0; k= l1) { - for (j=1; j idp) idlj -= idp; - war = wa[idlj - 2]; - wai = wa[idlj-1]; - for (ik=0; ik= l1) { - for (j=1; j= l1) { - for (k=0; k= l1) { - for (j=1; j= l1) { - for (k=0; k= l1) { - for (j=1; j= l1) { - for (j=1; j 5) { - wa[i1-1] = wa[i-1]; - wa[i1] = wa[i]; - } - } - l1 = l2; - } - } /* cffti1 */ - - -void npy_cffti(int n, Treal wsave[]) - { - int iw1, iw2; - if (n == 1) return; - iw1 = 2*n; - iw2 = iw1 + 2*n; - cffti1(n, wsave+iw1, (int*)(wsave+iw2)); - } /* npy_cffti */ - - /* ------------------------------------------------------------------- -rfftf1, rfftb1, npy_rfftf, npy_rfftb, rffti1, npy_rffti. Treal FFTs. ----------------------------------------------------------------------- */ - -static void rfftf1(int n, Treal c[], Treal ch[], const Treal wa[], const int ifac[MAXFAC+2]) - { - int i; - int k1, l1, l2, na, kh, nf, ip, iw, ix2, ix3, ix4, ido, idl1; - Treal *cinput, *coutput; - nf = ifac[1]; - na = 1; - l2 = n; - iw = n-1; - for (k1 = 1; k1 <= nf; ++k1) { - kh = nf - k1; - ip = ifac[kh + 2]; - l1 = l2 / ip; - ido = n / l2; - idl1 = ido*l1; - iw -= (ip - 1)*ido; - na = !na; - if (na) { - cinput = ch; - coutput = c; - } else { - cinput = c; - coutput = ch; - } - switch (ip) { - case 4: - ix2 = iw + ido; - ix3 = ix2 + ido; - radf4(ido, l1, cinput, coutput, &wa[iw], &wa[ix2], &wa[ix3]); - break; - case 2: - radf2(ido, l1, cinput, coutput, &wa[iw]); - break; - case 3: - ix2 = iw + ido; - radf3(ido, l1, cinput, coutput, &wa[iw], &wa[ix2]); - break; - case 5: - ix2 = iw + ido; - ix3 = ix2 + ido; - ix4 = ix3 + ido; - radf5(ido, l1, cinput, coutput, &wa[iw], &wa[ix2], &wa[ix3], &wa[ix4]); - break; - default: - if (ido == 1) - na = !na; - if (na == 0) { - radfg(ido, ip, l1, idl1, c, ch, &wa[iw]); - na = 1; - } else { - radfg(ido, ip, l1, idl1, ch, c, &wa[iw]); - na = 0; - } - } - l2 = l1; - } - if (na == 1) return; - for (i = 0; i < n; i++) c[i] = ch[i]; - } /* rfftf1 */ - - -static void rfftb1(int n, Treal c[], Treal ch[], const Treal wa[], const int ifac[MAXFAC+2]) - { - int i; - int k1, l1, l2, na, nf, ip, iw, ix2, ix3, ix4, ido, idl1; - Treal *cinput, *coutput; - nf = ifac[1]; - na = 0; - l1 = 1; - iw = 0; - for (k1=1; k1<=nf; k1++) { - ip = ifac[k1 + 1]; - l2 = ip*l1; - ido = n / l2; - idl1 = ido*l1; - if (na) { - cinput = ch; - coutput = c; - } else { - cinput = c; - coutput = ch; - } - switch (ip) { - case 4: - ix2 = iw + ido; - ix3 = ix2 + ido; - radb4(ido, l1, cinput, coutput, &wa[iw], &wa[ix2], &wa[ix3]); - na = !na; - break; - case 2: - radb2(ido, l1, cinput, coutput, &wa[iw]); - na = !na; - break; - case 3: - ix2 = iw + ido; - radb3(ido, l1, cinput, coutput, &wa[iw], &wa[ix2]); - na = !na; - break; - case 5: - ix2 = iw + ido; - ix3 = ix2 + ido; - ix4 = ix3 + ido; - radb5(ido, l1, cinput, coutput, &wa[iw], &wa[ix2], &wa[ix3], &wa[ix4]); - na = !na; - break; - default: - radbg(ido, ip, l1, idl1, cinput, coutput, &wa[iw]); - if (ido == 1) na = !na; - } - l1 = l2; - iw += (ip - 1)*ido; - } - if (na == 0) return; - for (i=0; i