Merge from branch bqresample

This commit is contained in:
Chris Cannam
2021-05-11 08:21:17 +01:00
26 changed files with 635 additions and 215 deletions

View File

@@ -11,3 +11,4 @@ efbc861f9b9460068c48a250232d343ffa7d5726 v1.7
d4911a276d96f6232a68c6b8448056d3946043b9 v1.8.1 d4911a276d96f6232a68c6b8448056d3946043b9 v1.8.1
fa6a54be7e6bf0c5adffd19ccec622481a8140a5 v1.8.2 fa6a54be7e6bf0c5adffd19ccec622481a8140a5 v1.8.2
37b18c17c042eafc39483ffb58837de844ba3289 v1.9 37b18c17c042eafc39483ffb58837de844ba3289 v1.9
7af7a76bbb1dc75f630555e22ca8f6ae9da79529 v1.9.1

View File

@@ -1,4 +1,24 @@
Changes in Rubber Band v1.9.1
* Switch build system from Makefiles and Visual Studio project to
Meson/Ninja for all platforms. There are still Makefiles and VS
projects included in otherbuilds/ for those who wish to use them to
build the static library directly
* Make various fixes to improve sound quality when pitch-shifting
dynamically in real-time (requires libsamplerate)
* Fix floating-point exception when a very very long stretch factor
is presented
* Move the two directories that together provide the .NET interface
(rubberband-sharp and rubberband-dll) into their own subdir (dotnet)
* Ensure the library builds and runs correctly on ARM Mac (Apple
Silicon, M1), and also on Windows using the Visual C++ Clang
front-end
The API is unchanged and the library is binary compatible with
version 1.7.
Changes in Rubber Band v1.9 Changes in Rubber Band v1.9
* Fix incorrect numbering of pitch speed/quality flags in the * Fix incorrect numbering of pitch speed/quality flags in the

View File

@@ -31,7 +31,7 @@ PROJECT_NAME = "Rubber Band Library"
# This could be handy for archiving the generated documentation or # This could be handy for archiving the generated documentation or
# if some version control system is used. # if some version control system is used.
PROJECT_NUMBER = 1.9.0 PROJECT_NUMBER = 1.9.1
# The OUTPUT_DIRECTORY tag is used to specify the (relative or absolute) # The OUTPUT_DIRECTORY tag is used to specify the (relative or absolute)
# base path where the generated documentation will be put. # base path where the generated documentation will be put.

View File

@@ -225,25 +225,28 @@ $ make -f otherbuilds/Makefile.linux
Ensure the Xcode command-line tools are installed, and ideally also Ensure the Xcode command-line tools are installed, and ideally also
install libsamplerate and libsndfile. install libsamplerate and libsndfile.
To build for the native architecture of the build machine: To build for the default architecture:
``` ```
$ meson build && ninja -C build $ meson build && ninja -C build
``` ```
To build for Intel (x86_64) regardless of the native architecture: Which architecture is the default may depend on the version of Meson
and/or the current shell. To force a particular architecture you can
use a Meson cross-file, as follows.
``` To build for Apple Silicon (arm64):
$ meson build --cross-file cross/macos-x86_64.txt && ninja -C build
```
To build for Apple Silicon (arm64) regardless of the native
architecture:
``` ```
$ meson build --cross-file cross/macos-arm64.txt && ninja -C build $ meson build --cross-file cross/macos-arm64.txt && ninja -C build
``` ```
To build for Intel (x86_64):
```
$ meson build --cross-file cross/macos-x86_64.txt && ninja -C build
```
You can build a universal binary library for both architectures like You can build a universal binary library for both architectures like
this: this:
@@ -317,6 +320,13 @@ The Rubber Band code is compatible with both the traditional Visual
C++ compiler (`cl`) and the Clang front-end (`clang`), and the build C++ compiler (`cl`) and the Clang front-end (`clang`), and the build
system will use whichever appears (first) in your path. system will use whichever appears (first) in your path.
To build against a specific Visual C++ runtime, use the built-in Meson
option `b_vscrt`:
```
> meson build -Db_vscrt=mt
```
See "FFT and resampler selection" below for further build options. See "FFT and resampler selection" below for further build options.
Alternatively, if you only need the static library and prefer a Visual Alternatively, if you only need the static library and prefer a Visual

View File

@@ -11,6 +11,7 @@ needs_exe_wrapper = false
c = 'cc' c = 'cc'
cpp = 'c++' cpp = 'c++'
strip = 'strip' strip = 'strip'
pkgconfig = 'pkg-config'
[built-in options] [built-in options]
c_args = ['-arch', 'arm64'] c_args = ['-arch', 'arm64']

View File

@@ -11,6 +11,7 @@ needs_exe_wrapper = false
c = 'cc' c = 'cc'
cpp = 'c++' cpp = 'c++'
strip = 'strip' strip = 'strip'
pkgconfig = 'pkg-config'
[built-in options] [built-in options]
c_args = ['-arch', 'arm64', '-arch', 'x86_64'] c_args = ['-arch', 'arm64', '-arch', 'x86_64']

View File

@@ -11,6 +11,7 @@ needs_exe_wrapper = false
c = 'cc' c = 'cc'
cpp = 'c++' cpp = 'c++'
strip = 'strip' strip = 'strip'
pkgconfig = 'pkg-config'
[built-in options] [built-in options]
c_args = ['-arch', 'x86_64'] c_args = ['-arch', 'x86_64']

View File

@@ -799,11 +799,19 @@ int main(int argc, char **argv)
usleep(10000); usleep(10000);
} }
} }
}
delete[] fbuf;
for (size_t i = 0; i < channels; ++i) delete[] ibuf[i];
delete[] ibuf;
}
sf_close(sndfile); sf_close(sndfile);
sf_close(sndfileOut); sf_close(sndfileOut);
free(fileName);
free(fileNameOut);
if (!quiet) { if (!quiet) {
cerr << "in: " << countIn << ", out: " << countOut << ", ratio: " << float(countOut)/float(countIn) << ", ideal output: " << lrint(countIn * ratio) << ", error: " << abs(lrint(countIn * ratio) - int(countOut)) << endl; cerr << "in: " << countIn << ", out: " << countOut << ", ratio: " << float(countOut)/float(countIn) << ", ideal output: " << lrint(countIn * ratio) << ", error: " << abs(lrint(countIn * ratio) - int(countOut)) << endl;
@@ -826,7 +834,7 @@ int main(int argc, char **argv)
} }
RubberBand::Profiler::dump(); RubberBand::Profiler::dump();
return 0; return 0;
} }

View File

@@ -2,7 +2,7 @@
project( project(
'Rubber Band Library', 'Rubber Band Library',
'c', 'cpp', 'c', 'cpp',
version: '1.9.0', version: '1.9.1',
license: 'GPL-2.0-or-later', license: 'GPL-2.0-or-later',
default_options: [ default_options: [
# All Rubber Band code is actually C++98, but some compilers no # All Rubber Band code is actually C++98, but some compilers no
@@ -15,9 +15,9 @@ project(
meson_version: '>= 0.53.0' meson_version: '>= 0.53.0'
) )
rubberband_dynamic_library_version = '2.1.2' rubberband_dynamic_library_version = '2.1.3'
system = build_machine.system() system = host_machine.system()
architecture = host_machine.cpu_family() architecture = host_machine.cpu_family()
cpp = meson.get_compiler('cpp') cpp = meson.get_compiler('cpp')
@@ -90,12 +90,21 @@ general_include_dirs = [
# Scan for any dependencies we may use later; all are optional # Scan for any dependencies we may use later; all are optional
extra_include_args = []
foreach d: get_option('extra_include_dirs')
extra_include_args += [ '-I' + d ]
endforeach
fftw3_dep = dependency('fftw3', version: '>= 3.0.0', required: false) fftw3_dep = dependency('fftw3', version: '>= 3.0.0', required: false)
samplerate_dep = dependency('samplerate', version: '>= 0.1.8', required: false) samplerate_dep = dependency('samplerate', version: '>= 0.1.8', required: false)
sndfile_dep = dependency('sndfile', version: '>= 1.0.16', required: false) sndfile_dep = dependency('sndfile', version: '>= 1.0.16', required: false)
vamp_dep = dependency('vamp-sdk', version: '>= 2.9', required: false) vamp_dep = dependency('vamp-sdk', version: '>= 2.9', required: false)
thread_dep = dependency('threads') thread_dep = dependency('threads')
have_ladspa = cpp.has_header('ladspa.h') have_ladspa = cpp.has_header('ladspa.h', args: extra_include_args)
have_jni = cpp.has_header('jni.h', args: extra_include_args)
javac = find_program('javac', required: false)
jar = find_program('jar', required: false)
# Check FFT and resampler options and set up dependencies and paths # Check FFT and resampler options and set up dependencies and paths
@@ -131,11 +140,6 @@ if resampler == 'auto'
endif endif
endif endif
extra_include_args = []
foreach d: get_option('extra_include_dirs')
extra_include_args += [ '-I' + d ]
endforeach
if fft == 'kissfft' if fft == 'kissfft'
config_summary += { 'FFT': 'KissFFT' } config_summary += { 'FFT': 'KissFFT' }
message('For FFT: using KissFFT') message('For FFT: using KissFFT')
@@ -362,6 +366,7 @@ if system == 'windows'
rubberband_program_name = 'rubberband-program' rubberband_program_name = 'rubberband-program'
rubberband_ladspa_name = 'ladspa-rubberband' rubberband_ladspa_name = 'ladspa-rubberband'
rubberband_vamp_name = 'vamp-rubberband' rubberband_vamp_name = 'vamp-rubberband'
rubberband_jni_name = 'rubberband-jni'
# Meson likes libxxx.a even on Windows, and so might we for a # Meson likes libxxx.a even on Windows, and so might we for a
# new library, but not here # new library, but not here
platform_static_name_prefix = '' platform_static_name_prefix = ''
@@ -372,13 +377,14 @@ else
rubberband_program_name = 'rubberband' rubberband_program_name = 'rubberband'
rubberband_ladspa_name = 'ladspa-rubberband' rubberband_ladspa_name = 'ladspa-rubberband'
rubberband_vamp_name = 'vamp-rubberband' rubberband_vamp_name = 'vamp-rubberband'
rubberband_jni_name = 'rubberband-jni'
platform_static_name_prefix = 'lib' platform_static_name_prefix = 'lib'
platform_static_name_suffix = 'a' platform_static_name_suffix = 'a'
endif endif
# And the build targets: Static and dynamic libraries, command-line # And the build targets: Static and dynamic libraries, command-line
# utility, LADSPA plugin, Vamp plugin # utility, LADSPA plugin, Vamp plugin, JNI library
message('Will build Rubber Band Library static library') message('Will build Rubber Band Library static library')
target_summary += { 'Static library': [ true, 'Name: ' + rubberband_static_name ] } target_summary += { 'Static library': [ true, 'Name: ' + rubberband_static_name ] }
@@ -419,6 +425,55 @@ else
message('Not building Rubber Band Library dynamic library: no_shared option set') message('Not building Rubber Band Library dynamic library: no_shared option set')
endif endif
if have_jni and javac.found() and jar.found()
target_summary += { 'JNI library': [ true, 'Name: ' + rubberband_jni_name ] }
message('Will build Java Native Interface')
rubberband_jni = shared_library(
rubberband_jni_name,
jni_sources,
include_directories: general_include_dirs,
cpp_args: general_compile_args,
c_args: general_compile_args,
link_args: [
arch_flags,
feature_libraries,
],
dependencies: [
rubberband_static_dep,
general_dependencies,
],
# NB the JNI library is not versioned
install: true,
)
rubberband_class = custom_target(
'rubberband_class',
input: 'com/breakfastquay/rubberband/RubberBandStretcher.java',
output: 'RubberBandStretcher.class',
command: [ javac, '@INPUT@', '-d', '@OUTDIR@' ],
)
rubberband_jar = custom_target(
'rubberband_jar',
input: rubberband_class,
output: 'rubberband.jar',
command: [ jar, 'cvf', '@OUTPUT@', 'com/breakfastquay/rubberband/@INPUT@' ],
build_by_default: true,
)
else
target_summary += { 'JNI library': false }
if not have_jni
message('Not building Java Native Interface: jni.h header not found')
else
message('Not building Java Native Interface: Java compiler not found')
endif
endif
install_headers(
[ 'rubberband/RubberBandStretcher.h',
'rubberband/rubberband-c.h'
],
subdir: 'rubberband'
)
if have_ladspa if have_ladspa
target_summary += { 'LADSPA plugin': [ true, 'Name: ' + rubberband_ladspa_name ] } target_summary += { 'LADSPA plugin': [ true, 'Name: ' + rubberband_ladspa_name ] }
message('Will build LADSPA plugin') message('Will build LADSPA plugin')

33
otherbuilds/deploy/macos.sh Executable file
View File

@@ -0,0 +1,33 @@
#!/bin/bash
set -eu
if [ ! -f ../rba/deploy/macos/notarize.sh ]; then
echo "need notarize script in ../rba/deploy/macos"
fi
version=$(grep '^ *version:' meson.build | head -1 | sed "s/^.*'\([0-9][0-9.]*\)'.*$/\1/")
echo
echo "Packaging command-line utility for Mac for Rubber Band v$version..."
echo
rm -rf build
PKG_CONFIG_PATH=/usr/local/lib/pkgconfig/ meson build --cross-file ./cross/macos-universal.txt
ninja -C build
./build/rubberband -V
key="Developer ID Application: Particular Programs Ltd (73F996B92S)"
mkdir -p packages
( cd build
codesign -s "$key" -fv --options runtime rubberband
zipfile="rubberband-$version-gpl-executable-macos.zip"
rm -f "$zipfile"
ditto -c -k rubberband "$zipfile"
../../rba/deploy/macos/notarize.sh "$zipfile" com.breakfastquay.rubberband
)
package_dir="rubberband-$version-gpl-executable-macos"
rm -rf "$package_dir"
mkdir "$package_dir"
cp build/rubberband "$package_dir"
cp CHANGELOG README.md COPYING "$package_dir"
tar cvjf "$package_dir.tar.bz2" "$package_dir"
mv "$package_dir.tar.bz2" packages/
rm -rf "$package_dir"
echo
echo "Done, package is in packages/$package_dir.tar.bz2"

59
otherbuilds/deploy/source.sh Executable file
View File

@@ -0,0 +1,59 @@
#!/bin/bash
set -eu
version=$(grep '^ *version:' meson.build | head -1 | sed "s/^.*'\([0-9][0-9.]*\)'.*$/\1/")
check() {
text="$1"
echo -n "$text [yN] "
read yn
case "$yn" in [Yy]) ;; *) echo "Exiting"; exit 3 ;; esac
}
echo
echo "Preparing to make source package for Rubber Band v$version..."
echo
grep '^ *version:' meson.build | head -1
check "Is the above version number (in meson.build) correct?"
echo
echo "The dynamic library version should have a point increment for each"
echo "release, a minor increment for backward-compatible ABI changes, and"
echo "a major increment for incompatible ABI changes."
echo
grep 'rubberband_dynamic_library_version' meson.build | head -1
check "Is the above library version (from meson.build) correct?"
echo
echo "The API major version should increment for incompatible API changes,"
echo "and the minor version should increment for backward-compatible API"
echo "changes."
echo
grep 'RUBBERBAND.*VERSION' rubberband/RubberBandStretcher.h
check "Are the above version and API versions (from the C++ header) correct?"
echo
echo "The C header should contain the same versions as the C++ header."
echo
grep 'RUBBERBAND.*VERSION' rubberband/rubberband-c.h
check "Are the above version and API versions (from the C header) correct?"
echo
grep '^PROJECT_NUMBER' Doxyfile
check "Is the above version (from Doxyfile) correct?"
echo
echo "Going ahead..."
mkdir -p packages
output="packages/rubberband-$version.tar.bz2"
hg archive --exclude otherbuilds/deploy "$output"
echo "Checking that the package compiles..."
tmpdir=$(mktemp -d)
cleanup() {
rm -rf "$tmpdir"
}
trap cleanup 0
prevdir=$(pwd)
( cd "$tmpdir"
tar xvf "$prevdir/$output"
cd "rubberband-$version"
meson build
ninja -C build
)
echo
echo "Checked, package is in $output"

View File

@@ -0,0 +1,54 @@
echo on
set STARTPWD=%CD%
set ORIGINALPATH=%PATH%
set PATH=C:\Program Files (x86)\Windows Kits\10\bin\x64;%PATH%
set vcvarsall="C:\Program Files (x86)\Microsoft Visual Studio\2019\Community\VC\Auxiliary\Build\vcvarsall.bat"
if not exist %vcvarsall% (
@ echo "Could not find MSVC vars batch file"
@ exit /b 2
)
set VERSION=%1
shift
if "%VERSION%" == "" (
@echo "Usage: win.bat <version>"
exit /b 1
)
@echo Building version %VERSION%
call %vcvarsall% amd64
if errorlevel 1 exit /b %errorlevel%
del /q /s build
meson build --buildtype release "-Dextra_include_dirs=C:\Program Files\libsndfile\include" "-Dextra_lib_dirs=C:\Program Files\libsndfile\lib" "-Db_vscrt=mt"
if errorlevel 1 exit /b %errorlevel%
ninja -C build
if errorlevel 1 exit /b %errorlevel%
cd build
ren rubberband-program.exe rubberband.exe
set NAME=Christopher Cannam
signtool sign /v /n "%NAME%" /t http://time.certum.pl /fd sha1 /a rubberband.exe
if errorlevel 1 exit /b %errorlevel%
cd ..
set DIR=rubberband-%VERSION%-gpl-executable-windows
del /q /s %DIR%
mkdir %DIR%
copy build\rubberband.exe %DIR%
copy "c:\Program Files\libsndfile\bin\sndfile.dll" %DIR%
copy COPYING %DIR%
copy README.md %DIR%
copy CHANGELOG %DIR%
set PATH=%ORIGINALPATH%
cd %STARTPWD%
@echo Done, now test and zip the directory %DIR%

View File

@@ -24,7 +24,7 @@
#ifndef RUBBERBAND_STRETCHER_H #ifndef RUBBERBAND_STRETCHER_H
#define RUBBERBAND_STRETCHER_H #define RUBBERBAND_STRETCHER_H
#define RUBBERBAND_VERSION "1.9.0" #define RUBBERBAND_VERSION "1.9.1"
#define RUBBERBAND_API_MAJOR_VERSION 2 #define RUBBERBAND_API_MAJOR_VERSION 2
#define RUBBERBAND_API_MINOR_VERSION 6 #define RUBBERBAND_API_MINOR_VERSION 6
@@ -186,7 +186,7 @@ public:
* one processing thread per audio channel in offline mode if * one processing thread per audio channel in offline mode if
* the stretcher is able to determine that more than one CPU is * the stretcher is able to determine that more than one CPU is
* available, and one thread only in realtime mode. This is the * available, and one thread only in realtime mode. This is the
* defafult. * default.
* *
* \li \c OptionThreadingNever - Never use more than one thread. * \li \c OptionThreadingNever - Never use more than one thread.
* *
@@ -464,7 +464,7 @@ public:
/** /**
* Change an OptionPitch configuration setting. This may be * Change an OptionPitch configuration setting. This may be
* called at any time in RealTime mode. It may not be called in * called at any time in RealTime mode. It may not be called in
* Offline mode (for which the transients option is fixed on * Offline mode (for which the pitch option is fixed on
* construction). * construction).
*/ */
void setPitchOption(Options options); void setPitchOption(Options options);

View File

@@ -28,7 +28,7 @@
extern "C" { extern "C" {
#endif #endif
#define RUBBERBAND_VERSION "1.9.0" #define RUBBERBAND_VERSION "1.9.1"
#define RUBBERBAND_API_MAJOR_VERSION 2 #define RUBBERBAND_API_MAJOR_VERSION 2
#define RUBBERBAND_API_MINOR_VERSION 6 #define RUBBERBAND_API_MINOR_VERSION 6

View File

@@ -58,14 +58,14 @@ HighFrequencyAudioCurve::processFloat(const float *R__ mag, int)
double double
HighFrequencyAudioCurve::processDouble(const double *R__ mag, int) HighFrequencyAudioCurve::processDouble(const double *R__ mag, int)
{ {
float result = 0.0; double result = 0.0;
const int sz = m_lastPerceivedBin; const int sz = m_lastPerceivedBin;
for (int n = 0; n <= sz; ++n) { for (int n = 0; n <= sz; ++n) {
result = result + mag[n] * n; result = result + mag[n] * n;
} }
return result; return result;
} }

View File

@@ -23,6 +23,8 @@
#include "Profiler.h" #include "Profiler.h"
#include "system/Thread.h"
#include <algorithm> #include <algorithm>
#include <set> #include <set>
#include <string> #include <string>
@@ -45,9 +47,13 @@ Profiler::m_profiles;
Profiler::WorstCallMap Profiler::WorstCallMap
Profiler::m_worstCalls; Profiler::m_worstCalls;
static Mutex profileMutex;
void void
Profiler::add(const char *id, float ms) Profiler::add(const char *id, float ms)
{ {
profileMutex.lock();
ProfileMap::iterator pmi = m_profiles.find(id); ProfileMap::iterator pmi = m_profiles.find(id);
if (pmi != m_profiles.end()) { if (pmi != m_profiles.end()) {
++pmi->second.first; ++pmi->second.first;
@@ -62,6 +68,8 @@ Profiler::add(const char *id, float ms)
} else { } else {
m_worstCalls[id] = ms; m_worstCalls[id] = ms;
} }
profileMutex.unlock();
} }
void void
@@ -74,6 +82,8 @@ Profiler::dump()
std::string std::string
Profiler::getReport() Profiler::getReport()
{ {
profileMutex.lock();
static const int buflen = 256; static const int buflen = 256;
char buffer[buflen]; char buffer[buflen];
std::string report; std::string report;
@@ -167,6 +177,8 @@ Profiler::getReport()
report += buffer; report += buffer;
} }
profileMutex.unlock();
return report; return report;
} }

View File

@@ -1851,7 +1851,7 @@ public:
dbuf[i] = realIn[i]; dbuf[i] = realIn[i];
} }
fftw_execute(m_dplanf); fftw_execute(m_dplanf);
v_convert(complexOut, (fft_double_type *)m_dpacked, sz + 2); v_convert(complexOut, (const fft_double_type *)m_dpacked, sz + 2);
} }
void forwardPolar(const double *R__ realIn, double *R__ magOut, double *R__ phaseOut) { void forwardPolar(const double *R__ realIn, double *R__ magOut, double *R__ phaseOut) {
@@ -1865,8 +1865,8 @@ public:
dbuf[i] = realIn[i]; dbuf[i] = realIn[i];
} }
fftw_execute(m_dplanf); fftw_execute(m_dplanf);
v_cartesian_interleaved_to_polar(magOut, phaseOut, v_cartesian_interleaved_to_polar
(double *)m_dpacked, m_size/2+1); (magOut, phaseOut, (const fft_double_type *)m_dpacked, m_size/2+1);
} }
void forwardMagnitude(const double *R__ realIn, double *R__ magOut) { void forwardMagnitude(const double *R__ realIn, double *R__ magOut) {
@@ -1880,11 +1880,8 @@ public:
dbuf[i] = realIn[i]; dbuf[i] = realIn[i];
} }
fftw_execute(m_dplanf); fftw_execute(m_dplanf);
const int hs = m_size/2; v_cartesian_interleaved_to_magnitudes
for (int i = 0; i <= hs; ++i) { (magOut, (const fft_double_type *)m_dpacked, m_size/2+1);
magOut[i] = sqrt(m_dpacked[i][0] * m_dpacked[i][0] +
m_dpacked[i][1] * m_dpacked[i][1]);
}
} }
void forward(const float *R__ realIn, float *R__ realOut, float *R__ imagOut) { void forward(const float *R__ realIn, float *R__ realOut, float *R__ imagOut) {
@@ -1912,7 +1909,7 @@ public:
fbuf[i] = realIn[i]; fbuf[i] = realIn[i];
} }
fftwf_execute(m_fplanf); fftwf_execute(m_fplanf);
v_convert(complexOut, (fft_float_type *)m_fpacked, sz + 2); v_convert(complexOut, (const fft_float_type *)m_fpacked, sz + 2);
} }
void forwardPolar(const float *R__ realIn, float *R__ magOut, float *R__ phaseOut) { void forwardPolar(const float *R__ realIn, float *R__ magOut, float *R__ phaseOut) {
@@ -1926,8 +1923,8 @@ public:
fbuf[i] = realIn[i]; fbuf[i] = realIn[i];
} }
fftwf_execute(m_fplanf); fftwf_execute(m_fplanf);
v_cartesian_interleaved_to_polar(magOut, phaseOut, v_cartesian_interleaved_to_polar
(float *)m_fpacked, m_size/2+1); (magOut, phaseOut, (fft_float_type *)m_fpacked, m_size/2+1);
} }
void forwardMagnitude(const float *R__ realIn, float *R__ magOut) { void forwardMagnitude(const float *R__ realIn, float *R__ magOut) {
@@ -1941,11 +1938,8 @@ public:
fbuf[i] = realIn[i]; fbuf[i] = realIn[i];
} }
fftwf_execute(m_fplanf); fftwf_execute(m_fplanf);
const int hs = m_size/2; v_cartesian_interleaved_to_magnitudes
for (int i = 0; i <= hs; ++i) { (magOut, (const fft_float_type *)m_fpacked, m_size/2+1);
magOut[i] = sqrtf(m_fpacked[i][0] * m_fpacked[i][0] +
m_fpacked[i][1] * m_fpacked[i][1]);
}
} }
void inverse(const double *R__ realIn, const double *R__ imagIn, double *R__ realOut) { void inverse(const double *R__ realIn, const double *R__ imagIn, double *R__ realOut) {
@@ -1964,7 +1958,7 @@ public:
void inverseInterleaved(const double *R__ complexIn, double *R__ realOut) { void inverseInterleaved(const double *R__ complexIn, double *R__ realOut) {
if (!m_dplanf) initDouble(); if (!m_dplanf) initDouble();
v_convert((double *)m_dpacked, complexIn, m_size + 2); v_convert((fft_double_type *)m_dpacked, complexIn, m_size + 2);
fftw_execute(m_dplani); fftw_execute(m_dplani);
const int sz = m_size; const int sz = m_size;
fft_double_type *const R__ dbuf = m_dbuf; fft_double_type *const R__ dbuf = m_dbuf;
@@ -1978,14 +1972,8 @@ public:
void inversePolar(const double *R__ magIn, const double *R__ phaseIn, double *R__ realOut) { void inversePolar(const double *R__ magIn, const double *R__ phaseIn, double *R__ realOut) {
if (!m_dplanf) initDouble(); if (!m_dplanf) initDouble();
const int hs = m_size/2; v_polar_to_cartesian_interleaved
fftw_complex *const R__ dpacked = m_dpacked; ((fft_double_type *)m_dpacked, magIn, phaseIn, m_size/2+1);
for (int i = 0; i <= hs; ++i) {
dpacked[i][0] = magIn[i] * cos(phaseIn[i]);
}
for (int i = 0; i <= hs; ++i) {
dpacked[i][1] = magIn[i] * sin(phaseIn[i]);
}
fftw_execute(m_dplani); fftw_execute(m_dplani);
const int sz = m_size; const int sz = m_size;
fft_double_type *const R__ dbuf = m_dbuf; fft_double_type *const R__ dbuf = m_dbuf;
@@ -2034,7 +2022,7 @@ public:
void inverseInterleaved(const float *R__ complexIn, float *R__ realOut) { void inverseInterleaved(const float *R__ complexIn, float *R__ realOut) {
if (!m_fplanf) initFloat(); if (!m_fplanf) initFloat();
v_copy((float *)m_fpacked, complexIn, m_size + 2); v_convert((fft_float_type *)m_fpacked, complexIn, m_size + 2);
fftwf_execute(m_fplani); fftwf_execute(m_fplani);
const int sz = m_size; const int sz = m_size;
fft_float_type *const R__ fbuf = m_fbuf; fft_float_type *const R__ fbuf = m_fbuf;
@@ -2048,14 +2036,8 @@ public:
void inversePolar(const float *R__ magIn, const float *R__ phaseIn, float *R__ realOut) { void inversePolar(const float *R__ magIn, const float *R__ phaseIn, float *R__ realOut) {
if (!m_fplanf) initFloat(); if (!m_fplanf) initFloat();
const int hs = m_size/2; v_polar_to_cartesian_interleaved
fftwf_complex *const R__ fpacked = m_fpacked; ((fft_float_type *)m_fpacked, magIn, phaseIn, m_size/2+1);
for (int i = 0; i <= hs; ++i) {
fpacked[i][0] = magIn[i] * cosf(phaseIn[i]);
}
for (int i = 0; i <= hs; ++i) {
fpacked[i][1] = magIn[i] * sinf(phaseIn[i]);
}
fftwf_execute(m_fplani); fftwf_execute(m_fplani);
const int sz = m_size; const int sz = m_size;
fft_float_type *const R__ fbuf = m_fbuf; fft_float_type *const R__ fbuf = m_fbuf;
@@ -3258,6 +3240,7 @@ FFT::~FFT()
void void
FFT::forward(const double *R__ realIn, double *R__ realOut, double *R__ imagOut) FFT::forward(const double *R__ realIn, double *R__ realOut, double *R__ imagOut)
{ {
Profiler profiler("FFT::forward");
CHECK_NOT_NULL(realIn); CHECK_NOT_NULL(realIn);
CHECK_NOT_NULL(realOut); CHECK_NOT_NULL(realOut);
CHECK_NOT_NULL(imagOut); CHECK_NOT_NULL(imagOut);
@@ -3267,6 +3250,7 @@ FFT::forward(const double *R__ realIn, double *R__ realOut, double *R__ imagOut)
void void
FFT::forwardInterleaved(const double *R__ realIn, double *R__ complexOut) FFT::forwardInterleaved(const double *R__ realIn, double *R__ complexOut)
{ {
Profiler profiler("FFT::forwardInterleaved");
CHECK_NOT_NULL(realIn); CHECK_NOT_NULL(realIn);
CHECK_NOT_NULL(complexOut); CHECK_NOT_NULL(complexOut);
d->forwardInterleaved(realIn, complexOut); d->forwardInterleaved(realIn, complexOut);
@@ -3275,6 +3259,7 @@ FFT::forwardInterleaved(const double *R__ realIn, double *R__ complexOut)
void void
FFT::forwardPolar(const double *R__ realIn, double *R__ magOut, double *R__ phaseOut) FFT::forwardPolar(const double *R__ realIn, double *R__ magOut, double *R__ phaseOut)
{ {
Profiler profiler("FFT::forwardPolar");
CHECK_NOT_NULL(realIn); CHECK_NOT_NULL(realIn);
CHECK_NOT_NULL(magOut); CHECK_NOT_NULL(magOut);
CHECK_NOT_NULL(phaseOut); CHECK_NOT_NULL(phaseOut);
@@ -3284,6 +3269,7 @@ FFT::forwardPolar(const double *R__ realIn, double *R__ magOut, double *R__ phas
void void
FFT::forwardMagnitude(const double *R__ realIn, double *R__ magOut) FFT::forwardMagnitude(const double *R__ realIn, double *R__ magOut)
{ {
Profiler profiler("FFT::forwardMagnitude");
CHECK_NOT_NULL(realIn); CHECK_NOT_NULL(realIn);
CHECK_NOT_NULL(magOut); CHECK_NOT_NULL(magOut);
d->forwardMagnitude(realIn, magOut); d->forwardMagnitude(realIn, magOut);
@@ -3292,6 +3278,7 @@ FFT::forwardMagnitude(const double *R__ realIn, double *R__ magOut)
void void
FFT::forward(const float *R__ realIn, float *R__ realOut, float *R__ imagOut) FFT::forward(const float *R__ realIn, float *R__ realOut, float *R__ imagOut)
{ {
Profiler profiler("FFT::forward[float]");
CHECK_NOT_NULL(realIn); CHECK_NOT_NULL(realIn);
CHECK_NOT_NULL(realOut); CHECK_NOT_NULL(realOut);
CHECK_NOT_NULL(imagOut); CHECK_NOT_NULL(imagOut);
@@ -3301,6 +3288,7 @@ FFT::forward(const float *R__ realIn, float *R__ realOut, float *R__ imagOut)
void void
FFT::forwardInterleaved(const float *R__ realIn, float *R__ complexOut) FFT::forwardInterleaved(const float *R__ realIn, float *R__ complexOut)
{ {
Profiler profiler("FFT::forwardInterleaved[float]");
CHECK_NOT_NULL(realIn); CHECK_NOT_NULL(realIn);
CHECK_NOT_NULL(complexOut); CHECK_NOT_NULL(complexOut);
d->forwardInterleaved(realIn, complexOut); d->forwardInterleaved(realIn, complexOut);
@@ -3309,6 +3297,7 @@ FFT::forwardInterleaved(const float *R__ realIn, float *R__ complexOut)
void void
FFT::forwardPolar(const float *R__ realIn, float *R__ magOut, float *R__ phaseOut) FFT::forwardPolar(const float *R__ realIn, float *R__ magOut, float *R__ phaseOut)
{ {
Profiler profiler("FFT::forwardPolar[float]");
CHECK_NOT_NULL(realIn); CHECK_NOT_NULL(realIn);
CHECK_NOT_NULL(magOut); CHECK_NOT_NULL(magOut);
CHECK_NOT_NULL(phaseOut); CHECK_NOT_NULL(phaseOut);
@@ -3318,6 +3307,7 @@ FFT::forwardPolar(const float *R__ realIn, float *R__ magOut, float *R__ phaseOu
void void
FFT::forwardMagnitude(const float *R__ realIn, float *R__ magOut) FFT::forwardMagnitude(const float *R__ realIn, float *R__ magOut)
{ {
Profiler profiler("FFT::forwardMagnitude[float]");
CHECK_NOT_NULL(realIn); CHECK_NOT_NULL(realIn);
CHECK_NOT_NULL(magOut); CHECK_NOT_NULL(magOut);
d->forwardMagnitude(realIn, magOut); d->forwardMagnitude(realIn, magOut);
@@ -3326,6 +3316,7 @@ FFT::forwardMagnitude(const float *R__ realIn, float *R__ magOut)
void void
FFT::inverse(const double *R__ realIn, const double *R__ imagIn, double *R__ realOut) FFT::inverse(const double *R__ realIn, const double *R__ imagIn, double *R__ realOut)
{ {
Profiler profiler("FFT::inverse");
CHECK_NOT_NULL(realIn); CHECK_NOT_NULL(realIn);
CHECK_NOT_NULL(imagIn); CHECK_NOT_NULL(imagIn);
CHECK_NOT_NULL(realOut); CHECK_NOT_NULL(realOut);
@@ -3335,6 +3326,7 @@ FFT::inverse(const double *R__ realIn, const double *R__ imagIn, double *R__ rea
void void
FFT::inverseInterleaved(const double *R__ complexIn, double *R__ realOut) FFT::inverseInterleaved(const double *R__ complexIn, double *R__ realOut)
{ {
Profiler profiler("FFT::inverseInterleaved");
CHECK_NOT_NULL(complexIn); CHECK_NOT_NULL(complexIn);
CHECK_NOT_NULL(realOut); CHECK_NOT_NULL(realOut);
d->inverseInterleaved(complexIn, realOut); d->inverseInterleaved(complexIn, realOut);
@@ -3343,6 +3335,7 @@ FFT::inverseInterleaved(const double *R__ complexIn, double *R__ realOut)
void void
FFT::inversePolar(const double *R__ magIn, const double *R__ phaseIn, double *R__ realOut) FFT::inversePolar(const double *R__ magIn, const double *R__ phaseIn, double *R__ realOut)
{ {
Profiler profiler("FFT::inversePolar");
CHECK_NOT_NULL(magIn); CHECK_NOT_NULL(magIn);
CHECK_NOT_NULL(phaseIn); CHECK_NOT_NULL(phaseIn);
CHECK_NOT_NULL(realOut); CHECK_NOT_NULL(realOut);
@@ -3352,6 +3345,7 @@ FFT::inversePolar(const double *R__ magIn, const double *R__ phaseIn, double *R_
void void
FFT::inverseCepstral(const double *R__ magIn, double *R__ cepOut) FFT::inverseCepstral(const double *R__ magIn, double *R__ cepOut)
{ {
Profiler profiler("FFT::inverseCepstral");
CHECK_NOT_NULL(magIn); CHECK_NOT_NULL(magIn);
CHECK_NOT_NULL(cepOut); CHECK_NOT_NULL(cepOut);
d->inverseCepstral(magIn, cepOut); d->inverseCepstral(magIn, cepOut);
@@ -3360,6 +3354,7 @@ FFT::inverseCepstral(const double *R__ magIn, double *R__ cepOut)
void void
FFT::inverse(const float *R__ realIn, const float *R__ imagIn, float *R__ realOut) FFT::inverse(const float *R__ realIn, const float *R__ imagIn, float *R__ realOut)
{ {
Profiler profiler("FFT::inverse[float]");
CHECK_NOT_NULL(realIn); CHECK_NOT_NULL(realIn);
CHECK_NOT_NULL(imagIn); CHECK_NOT_NULL(imagIn);
CHECK_NOT_NULL(realOut); CHECK_NOT_NULL(realOut);
@@ -3369,6 +3364,7 @@ FFT::inverse(const float *R__ realIn, const float *R__ imagIn, float *R__ realOu
void void
FFT::inverseInterleaved(const float *R__ complexIn, float *R__ realOut) FFT::inverseInterleaved(const float *R__ complexIn, float *R__ realOut)
{ {
Profiler profiler("FFT::inverseInterleaved[float]");
CHECK_NOT_NULL(complexIn); CHECK_NOT_NULL(complexIn);
CHECK_NOT_NULL(realOut); CHECK_NOT_NULL(realOut);
d->inverseInterleaved(complexIn, realOut); d->inverseInterleaved(complexIn, realOut);
@@ -3377,6 +3373,7 @@ FFT::inverseInterleaved(const float *R__ complexIn, float *R__ realOut)
void void
FFT::inversePolar(const float *R__ magIn, const float *R__ phaseIn, float *R__ realOut) FFT::inversePolar(const float *R__ magIn, const float *R__ phaseIn, float *R__ realOut)
{ {
Profiler profiler("FFT::inversePolar[float]");
CHECK_NOT_NULL(magIn); CHECK_NOT_NULL(magIn);
CHECK_NOT_NULL(phaseIn); CHECK_NOT_NULL(phaseIn);
CHECK_NOT_NULL(realOut); CHECK_NOT_NULL(realOut);
@@ -3386,6 +3383,7 @@ FFT::inversePolar(const float *R__ magIn, const float *R__ phaseIn, float *R__ r
void void
FFT::inverseCepstral(const float *R__ magIn, float *R__ cepOut) FFT::inverseCepstral(const float *R__ magIn, float *R__ cepOut)
{ {
Profiler profiler("FFT::inverseCepstral[float]");
CHECK_NOT_NULL(magIn); CHECK_NOT_NULL(magIn);
CHECK_NOT_NULL(cepOut); CHECK_NOT_NULL(cepOut);
d->inverseCepstral(magIn, cepOut); d->inverseCepstral(magIn, cepOut);

View File

@@ -329,7 +329,6 @@ Java_com_breakfastquay_rubberband_RubberBandStretcher_study(JNIEnv *env, jobject
int channels = env->GetArrayLength(data); int channels = env->GetArrayLength(data);
float **arr = allocate<float *>(channels); float **arr = allocate<float *>(channels);
float **input = allocate<float *>(channels); float **input = allocate<float *>(channels);
int samples = 0;
for (int c = 0; c < channels; ++c) { for (int c = 0; c < channels; ++c) {
jfloatArray cdata = (jfloatArray)env->GetObjectArrayElement(data, c); jfloatArray cdata = (jfloatArray)env->GetObjectArrayElement(data, c);
arr[c] = env->GetFloatArrayElements(cdata, 0); arr[c] = env->GetFloatArrayElements(cdata, 0);
@@ -350,7 +349,6 @@ Java_com_breakfastquay_rubberband_RubberBandStretcher_process(JNIEnv *env, jobje
int channels = env->GetArrayLength(data); int channels = env->GetArrayLength(data);
float **arr = allocate<float *>(channels); float **arr = allocate<float *>(channels);
float **input = allocate<float *>(channels); float **input = allocate<float *>(channels);
int samples = 0;
for (int c = 0; c < channels; ++c) { for (int c = 0; c < channels; ++c) {
jfloatArray cdata = (jfloatArray)env->GetObjectArrayElement(data, c); jfloatArray cdata = (jfloatArray)env->GetObjectArrayElement(data, c);
arr[c] = env->GetFloatArrayElements(cdata, 0); arr[c] = env->GetFloatArrayElements(cdata, 0);
@@ -383,7 +381,7 @@ Java_com_breakfastquay_rubberband_RubberBandStretcher_retrieve(JNIEnv *env, jobj
float **outbuf = allocate_channels<float>(channels, n); float **outbuf = allocate_channels<float>(channels, n);
size_t retrieved = stretcher->retrieve(outbuf, n); size_t retrieved = stretcher->retrieve(outbuf, n);
for (int c = 0; c < channels; ++c) { for (size_t c = 0; c < channels; ++c) {
jfloatArray cdata = (jfloatArray)env->GetObjectArrayElement(output, c); jfloatArray cdata = (jfloatArray)env->GetObjectArrayElement(output, c);
env->SetFloatArrayRegion(cdata, offset, retrieved, outbuf[c]); env->SetFloatArrayRegion(cdata, offset, retrieved, outbuf[c]);
} }

View File

@@ -1,11 +1,11 @@
Copyright (c) 2003-2004 Mark Borgerding Copyright (c) 2003-2010 Mark Borgerding . All rights reserved.
All rights reserved. KISS FFT is provided under:
Redistribution and use in source and binary forms, with or without modification, are permitted provided that the following conditions are met: SPDX-License-Identifier: BSD-3-Clause
* Redistributions of source code must retain the above copyright notice, this list of conditions and the following disclaimer. Being under the terms of the BSD 3-clause "New" or "Revised" License,
* Redistributions in binary form must reproduce the above copyright notice, this list of conditions and the following disclaimer in the documentation and/or other materials provided with the distribution. according with:
* Neither the author nor the names of any contributors may be used to endorse or promote products derived from this software without specific prior written permission.
LICENSES/BSD-3-Clause
THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.

View File

@@ -1,29 +1,25 @@
/* /*
Copyright (c) 2003-2004, Mark Borgerding * Copyright (c) 2003-2010, Mark Borgerding. All rights reserved.
* This file is part of KISS FFT - https://github.com/mborgerding/kissfft
All rights reserved. *
* SPDX-License-Identifier: BSD-3-Clause
Redistribution and use in source and binary forms, with or without modification, are permitted provided that the following conditions are met: * See COPYING file for more information.
*/
* Redistributions of source code must retain the above copyright notice, this list of conditions and the following disclaimer.
* Redistributions in binary form must reproduce the above copyright notice, this list of conditions and the following disclaimer in the documentation and/or other materials provided with the distribution.
* Neither the author nor the names of any contributors may be used to endorse or promote products derived from this software without specific prior written permission.
THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
*/
#ifndef KISS_FFT_GUTS_H
#define KISS_FFT_GUTS_H
/* kiss_fft.h /* kiss_fft.h
defines kiss_fft_scalar as either short or a float type defines kiss_fft_scalar as either short or a float type
and defines and defines
typedef struct { kiss_fft_scalar r; kiss_fft_scalar i; }kiss_fft_cpx; */ typedef struct { kiss_fft_scalar r; kiss_fft_scalar i; }kiss_fft_cpx; */
#ifndef _kiss_fft_guts_h
#define _kiss_fft_guts_h
#include "kiss_fft.h" #include "kiss_fft.h"
#include "kiss_fft_log.h"
#include <limits.h> #include <limits.h>
#define MAXFACTORS 32 #define MAXFACTORS 32
/* e.g. an fft of length 128 has 4 factors /* e.g. an fft of length 128 has 4 factors
as far as kissfft is concerned as far as kissfft is concerned
4*4*4*2 4*4*4*2
*/ */
@@ -45,22 +41,23 @@ struct kiss_fft_state{
C_ADDTO( res , a) : res += a C_ADDTO( res , a) : res += a
* */ * */
#ifdef FIXED_POINT #ifdef FIXED_POINT
#include <stdint.h>
#if (FIXED_POINT==32) #if (FIXED_POINT==32)
# define FRACBITS 31 # define FRACBITS 31
# define SAMPPROD int64_t # define SAMPPROD int64_t
#define SAMP_MAX 2147483647 #define SAMP_MAX INT32_MAX
#define SAMP_MIN INT32_MIN
#else #else
# define FRACBITS 15 # define FRACBITS 15
# define SAMPPROD int32_t # define SAMPPROD int32_t
#define SAMP_MAX 32767 #define SAMP_MAX INT16_MAX
#define SAMP_MIN INT16_MIN
#endif #endif
#define SAMP_MIN -SAMP_MAX
#if defined(CHECK_OVERFLOW) #if defined(CHECK_OVERFLOW)
# define CHECK_OVERFLOW_OP(a,op,b) \ # define CHECK_OVERFLOW_OP(a,op,b) \
if ( (SAMPPROD)(a) op (SAMPPROD)(b) > SAMP_MAX || (SAMPPROD)(a) op (SAMPPROD)(b) < SAMP_MIN ) { \ if ( (SAMPPROD)(a) op (SAMPPROD)(b) > SAMP_MAX || (SAMPPROD)(a) op (SAMPPROD)(b) < SAMP_MIN ) { \
fprintf(stderr,"WARNING:overflow @ " __FILE__ "(%d): (%d " #op" %d) = %ld\n",__LINE__,(a),(b),(SAMPPROD)(a) op (SAMPPROD)(b) ); } KISS_FFT_WARNING("overflow (%d " #op" %d) = %ld", (a),(b),(SAMPPROD)(a) op (SAMPPROD)(b)); }
#endif #endif
@@ -74,11 +71,11 @@ struct kiss_fft_state{
(m).i = sround( smul((a).r,(b).i) + smul((a).i,(b).r) ); }while(0) (m).i = sround( smul((a).r,(b).i) + smul((a).i,(b).r) ); }while(0)
# define DIVSCALAR(x,k) \ # define DIVSCALAR(x,k) \
(x) = sround( smul( x, SAMP_MAX/k ) ) (x) = sround( smul( x, SAMP_MAX/k ) )
# define C_FIXDIV(c,div) \ # define C_FIXDIV(c,div) \
do { DIVSCALAR( (c).r , div); \ do { DIVSCALAR( (c).r , div); \
DIVSCALAR( (c).i , div); }while (0) DIVSCALAR( (c).i , div); }while (0)
# define C_MULBYSCALAR( c, s ) \ # define C_MULBYSCALAR( c, s ) \
do{ (c).r = sround( smul( (c).r , s ) ) ;\ do{ (c).r = sround( smul( (c).r , s ) ) ;\
@@ -102,28 +99,28 @@ struct kiss_fft_state{
#define C_ADD( res, a,b)\ #define C_ADD( res, a,b)\
do { \ do { \
CHECK_OVERFLOW_OP((a).r,+,(b).r)\ CHECK_OVERFLOW_OP((a).r,+,(b).r)\
CHECK_OVERFLOW_OP((a).i,+,(b).i)\ CHECK_OVERFLOW_OP((a).i,+,(b).i)\
(res).r=(a).r+(b).r; (res).i=(a).i+(b).i; \ (res).r=(a).r+(b).r; (res).i=(a).i+(b).i; \
}while(0) }while(0)
#define C_SUB( res, a,b)\ #define C_SUB( res, a,b)\
do { \ do { \
CHECK_OVERFLOW_OP((a).r,-,(b).r)\ CHECK_OVERFLOW_OP((a).r,-,(b).r)\
CHECK_OVERFLOW_OP((a).i,-,(b).i)\ CHECK_OVERFLOW_OP((a).i,-,(b).i)\
(res).r=(a).r-(b).r; (res).i=(a).i-(b).i; \ (res).r=(a).r-(b).r; (res).i=(a).i-(b).i; \
}while(0) }while(0)
#define C_ADDTO( res , a)\ #define C_ADDTO( res , a)\
do { \ do { \
CHECK_OVERFLOW_OP((res).r,+,(a).r)\ CHECK_OVERFLOW_OP((res).r,+,(a).r)\
CHECK_OVERFLOW_OP((res).i,+,(a).i)\ CHECK_OVERFLOW_OP((res).i,+,(a).i)\
(res).r += (a).r; (res).i += (a).i;\ (res).r += (a).r; (res).i += (a).i;\
}while(0) }while(0)
#define C_SUBFROM( res , a)\ #define C_SUBFROM( res , a)\
do {\ do {\
CHECK_OVERFLOW_OP((res).r,-,(a).r)\ CHECK_OVERFLOW_OP((res).r,-,(a).r)\
CHECK_OVERFLOW_OP((res).i,-,(a).i)\ CHECK_OVERFLOW_OP((res).i,-,(a).i)\
(res).r -= (a).r; (res).i -= (a).i; \ (res).r -= (a).r; (res).i -= (a).i; \
}while(0) }while(0)
@@ -138,18 +135,33 @@ struct kiss_fft_state{
#else #else
# define KISS_FFT_COS(phase) (kiss_fft_scalar) cos(phase) # define KISS_FFT_COS(phase) (kiss_fft_scalar) cos(phase)
# define KISS_FFT_SIN(phase) (kiss_fft_scalar) sin(phase) # define KISS_FFT_SIN(phase) (kiss_fft_scalar) sin(phase)
# define HALF_OF(x) ((x)*.5) # define HALF_OF(x) ((x)*((kiss_fft_scalar).5))
#endif #endif
#define kf_cexp(x,phase) \ #define kf_cexp(x,phase) \
do{ \ do{ \
(x)->r = KISS_FFT_COS(phase);\ (x)->r = KISS_FFT_COS(phase);\
(x)->i = KISS_FFT_SIN(phase);\ (x)->i = KISS_FFT_SIN(phase);\
}while(0) }while(0)
/* a debugging function */ /* a debugging function */
#define pcpx(c)\ #define pcpx(c)\
fprintf(stderr,"%g + %gi\n",(double)((c)->r),(double)((c)->i) ) KISS_FFT_DEBUG("%g + %gi\n",(double)((c)->r),(double)((c)->i))
#ifdef KISS_FFT_USE_ALLOCA
// define this to allow use of alloca instead of malloc for temporary buffers
// Temporary buffers are used in two case:
// 1. FFT sizes that have "bad" factors. i.e. not 2,3 and 5
// 2. "in-place" FFTs. Notice the quotes, since kissfft does not really do an in-place transform.
#include <alloca.h>
#define KISS_FFT_TMP_ALLOC(nbytes) alloca(nbytes)
#define KISS_FFT_TMP_FREE(ptr)
#else
#define KISS_FFT_TMP_ALLOC(nbytes) KISS_FFT_MALLOC(nbytes)
#define KISS_FFT_TMP_FREE(ptr) KISS_FFT_FREE(ptr)
#endif #endif
#endif /* _kiss_fft_guts_h */

View File

@@ -1,16 +1,10 @@
/* /*
Copyright (c) 2003-2004, Mark Borgerding * Copyright (c) 2003-2010, Mark Borgerding. All rights reserved.
* This file is part of KISS FFT - https://github.com/mborgerding/kissfft
All rights reserved. *
* SPDX-License-Identifier: BSD-3-Clause
Redistribution and use in source and binary forms, with or without modification, are permitted provided that the following conditions are met: * See COPYING file for more information.
*/
* Redistributions of source code must retain the above copyright notice, this list of conditions and the following disclaimer.
* Redistributions in binary form must reproduce the above copyright notice, this list of conditions and the following disclaimer in the documentation and/or other materials provided with the distribution.
* Neither the author nor the names of any contributors may be used to endorse or promote products derived from this software without specific prior written permission.
THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
*/
#include "_kiss_fft_guts.h" #include "_kiss_fft_guts.h"
@@ -18,21 +12,6 @@ THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND
fixed or floating point complex numbers. It also delares the kf_ internal functions. fixed or floating point complex numbers. It also delares the kf_ internal functions.
*/ */
static kiss_fft_cpx *scratchbuf=NULL;
static size_t nscratchbuf=0;
static kiss_fft_cpx *tmpbuf=NULL;
static size_t ntmpbuf=0;
#define CHECKBUF(buf,nbuf,n) \
do { \
if ( nbuf < (size_t)(n) ) {\
free(buf); \
buf = (kiss_fft_cpx*)KISS_FFT_MALLOC(sizeof(kiss_fft_cpx)*(n)); \
nbuf = (size_t)(n); \
} \
}while(0)
static void kf_bfly2( static void kf_bfly2(
kiss_fft_cpx * Fout, kiss_fft_cpx * Fout,
const size_t fstride, const size_t fstride,
@@ -69,6 +48,7 @@ static void kf_bfly4(
const size_t m2=2*m; const size_t m2=2*m;
const size_t m3=3*m; const size_t m3=3*m;
tw3 = tw2 = tw1 = st->twiddles; tw3 = tw2 = tw1 = st->twiddles;
do { do {
@@ -222,29 +202,34 @@ static void kf_bfly_generic(
kiss_fft_cpx t; kiss_fft_cpx t;
int Norig = st->nfft; int Norig = st->nfft;
CHECKBUF(scratchbuf,nscratchbuf,p); kiss_fft_cpx * scratch = (kiss_fft_cpx*)KISS_FFT_TMP_ALLOC(sizeof(kiss_fft_cpx)*p);
if (scratch == NULL){
KISS_FFT_ERROR("Memory allocation failed.");
return;
}
for ( u=0; u<m; ++u ) { for ( u=0; u<m; ++u ) {
k=u; k=u;
for ( q1=0 ; q1<p ; ++q1 ) { for ( q1=0 ; q1<p ; ++q1 ) {
scratchbuf[q1] = Fout[ k ]; scratch[q1] = Fout[ k ];
C_FIXDIV(scratchbuf[q1],p); C_FIXDIV(scratch[q1],p);
k += m; k += m;
} }
k=u; k=u;
for ( q1=0 ; q1<p ; ++q1 ) { for ( q1=0 ; q1<p ; ++q1 ) {
int twidx=0; int twidx=0;
Fout[ k ] = scratchbuf[0]; Fout[ k ] = scratch[0];
for (q=1;q<p;++q ) { for (q=1;q<p;++q ) {
twidx += fstride * k; twidx += fstride * k;
if (twidx>=Norig) twidx-=Norig; if (twidx>=Norig) twidx-=Norig;
C_MUL(t,scratchbuf[q] , twiddles[twidx] ); C_MUL(t,scratch[q] , twiddles[twidx] );
C_ADDTO( Fout[ k ] ,t); C_ADDTO( Fout[ k ] ,t);
} }
k += m; k += m;
} }
} }
KISS_FFT_TMP_FREE(scratch);
} }
static static
@@ -262,6 +247,30 @@ void kf_work(
const int m=*factors++; /* stage's fft length/p */ const int m=*factors++; /* stage's fft length/p */
const kiss_fft_cpx * Fout_end = Fout + p*m; const kiss_fft_cpx * Fout_end = Fout + p*m;
#ifdef _OPENMP
// use openmp extensions at the
// top-level (not recursive)
if (fstride==1 && p<=5 && m!=1)
{
int k;
// execute the p different work units in different threads
# pragma omp parallel for
for (k=0;k<p;++k)
kf_work( Fout +k*m, f+ fstride*in_stride*k,fstride*p,in_stride,factors,st);
// all threads have joined by this point
switch (p) {
case 2: kf_bfly2(Fout,fstride,st,m); break;
case 3: kf_bfly3(Fout,fstride,st,m); break;
case 4: kf_bfly4(Fout,fstride,st,m); break;
case 5: kf_bfly5(Fout,fstride,st,m); break;
default: kf_bfly_generic(Fout,fstride,st,m,p); break;
}
return;
}
#endif
if (m==1) { if (m==1) {
do{ do{
*Fout = *f; *Fout = *f;
@@ -269,6 +278,10 @@ void kf_work(
}while(++Fout != Fout_end ); }while(++Fout != Fout_end );
}else{ }else{
do{ do{
// recursive call:
// DFT of size m*p performed by doing
// p instances of smaller DFTs of size m,
// each one takes a decimated version of the input
kf_work( Fout , f, fstride*p, in_stride, factors,st); kf_work( Fout , f, fstride*p, in_stride, factors,st);
f += fstride*in_stride; f += fstride*in_stride;
}while( (Fout += m) != Fout_end ); }while( (Fout += m) != Fout_end );
@@ -276,20 +289,21 @@ void kf_work(
Fout=Fout_beg; Fout=Fout_beg;
// recombine the p smaller DFTs
switch (p) { switch (p) {
case 2: kf_bfly2(Fout,fstride,st,m); break; case 2: kf_bfly2(Fout,fstride,st,m); break;
case 3: kf_bfly3(Fout,fstride,st,m); break; case 3: kf_bfly3(Fout,fstride,st,m); break;
case 4: kf_bfly4(Fout,fstride,st,m); break; case 4: kf_bfly4(Fout,fstride,st,m); break;
case 5: kf_bfly5(Fout,fstride,st,m); break; case 5: kf_bfly5(Fout,fstride,st,m); break;
default: kf_bfly_generic(Fout,fstride,st,m,p); break; default: kf_bfly_generic(Fout,fstride,st,m,p); break;
} }
} }
/* facbuf is populated by p1,m1,p2,m2, ... /* facbuf is populated by p1,m1,p2,m2, ...
where where
p[i] * m[i] = m[i-1] p[i] * m[i] = m[i-1]
m0 = n */ m0 = n */
static static
void kf_factor(int n,int * facbuf) void kf_factor(int n,int * facbuf)
{ {
int p=4; int p=4;
@@ -322,9 +336,11 @@ void kf_factor(int n,int * facbuf)
* */ * */
kiss_fft_cfg kiss_fft_alloc(int nfft,int inverse_fft,void * mem,size_t * lenmem ) kiss_fft_cfg kiss_fft_alloc(int nfft,int inverse_fft,void * mem,size_t * lenmem )
{ {
KISS_FFT_ALIGN_CHECK(mem)
kiss_fft_cfg st=NULL; kiss_fft_cfg st=NULL;
size_t memneeded = sizeof(struct kiss_fft_state) size_t memneeded = KISS_FFT_ALIGN_SIZE_UP(sizeof(struct kiss_fft_state)
+ sizeof(kiss_fft_cpx)*(nfft-1); /* twiddle factors*/ + sizeof(kiss_fft_cpx)*(nfft-1)); /* twiddle factors*/
if ( lenmem==NULL ) { if ( lenmem==NULL ) {
st = ( kiss_fft_cfg)KISS_FFT_MALLOC( memneeded ); st = ( kiss_fft_cfg)KISS_FFT_MALLOC( memneeded );
@@ -352,14 +368,27 @@ kiss_fft_cfg kiss_fft_alloc(int nfft,int inverse_fft,void * mem,size_t * lenmem
} }
void kiss_fft_stride(kiss_fft_cfg st,const kiss_fft_cpx *fin,kiss_fft_cpx *fout,int in_stride) void kiss_fft_stride(kiss_fft_cfg st,const kiss_fft_cpx *fin,kiss_fft_cpx *fout,int in_stride)
{ {
if (fin == fout) { if (fin == fout) {
CHECKBUF(tmpbuf,ntmpbuf,st->nfft); //NOTE: this is not really an in-place FFT algorithm.
//It just performs an out-of-place FFT into a temp buffer
if (fout == NULL){
KISS_FFT_ERROR("fout buffer NULL.");
return;
}
kiss_fft_cpx * tmpbuf = (kiss_fft_cpx*)KISS_FFT_TMP_ALLOC( sizeof(kiss_fft_cpx)*st->nfft);
if (tmpbuf == NULL){
KISS_FFT_ERROR("Memory allocation error.");
return;
}
kf_work(tmpbuf,fin,1,in_stride, st->factors,st); kf_work(tmpbuf,fin,1,in_stride, st->factors,st);
memcpy(fout,tmpbuf,sizeof(kiss_fft_cpx)*st->nfft); memcpy(fout,tmpbuf,sizeof(kiss_fft_cpx)*st->nfft);
KISS_FFT_TMP_FREE(tmpbuf);
}else{ }else{
kf_work( fout, fin, 1,in_stride, st->factors,st ); kf_work( fout, fin, 1,in_stride, st->factors,st );
} }
@@ -371,17 +400,9 @@ void kiss_fft(kiss_fft_cfg cfg,const kiss_fft_cpx *fin,kiss_fft_cpx *fout)
} }
/* not really necessary to call, but if someone is doing in-place ffts, they may want to free the
buffers from CHECKBUF
*/
void kiss_fft_cleanup(void) void kiss_fft_cleanup(void)
{ {
free(scratchbuf); // nothing needed any more
scratchbuf = NULL;
nscratchbuf=0;
free(tmpbuf);
tmpbuf=NULL;
ntmpbuf=0;
} }
int kiss_fft_next_fast_size(int n) int kiss_fft_next_fast_size(int n)

View File

@@ -1,12 +1,32 @@
/*
* Copyright (c) 2003-2010, Mark Borgerding. All rights reserved.
* This file is part of KISS FFT - https://github.com/mborgerding/kissfft
*
* SPDX-License-Identifier: BSD-3-Clause
* See COPYING file for more information.
*/
#ifndef KISS_FFT_H #ifndef KISS_FFT_H
#define KISS_FFT_H #define KISS_FFT_H
#include <stdlib.h> #include <stdlib.h>
#include <stdio.h> #include <stdio.h>
#include <math.h> #include <math.h>
#include <memory.h> #include <string.h>
#ifndef __APPLE__
#include <malloc.h> // Define KISS_FFT_SHARED macro to properly export symbols
#ifdef KISS_FFT_SHARED
# ifdef _WIN32
# ifdef KISS_FFT_BUILD
# define KISS_FFT_API __declspec(dllexport)
# else
# define KISS_FFT_API __declspec(dllimport)
# endif
# else
# define KISS_FFT_API __attribute__ ((visibility ("default")))
# endif
#else
# define KISS_FFT_API
#endif #endif
#ifdef __cplusplus #ifdef __cplusplus
@@ -26,17 +46,32 @@ extern "C" {
in the tools/ directory. in the tools/ directory.
*/ */
/* User may override KISS_FFT_MALLOC and/or KISS_FFT_FREE. */
#ifdef USE_SIMD #ifdef USE_SIMD
# include <xmmintrin.h> # include <xmmintrin.h>
# define kiss_fft_scalar __m128 # define kiss_fft_scalar __m128
#define KISS_FFT_MALLOC(nbytes) memalign(16,nbytes) # ifndef KISS_FFT_MALLOC
#else # define KISS_FFT_MALLOC(nbytes) _mm_malloc(nbytes,16)
#define KISS_FFT_MALLOC malloc # define KISS_FFT_ALIGN_CHECK(ptr)
#endif # define KISS_FFT_ALIGN_SIZE_UP(size) ((size + 15UL) & ~0xFUL)
# endif
# ifndef KISS_FFT_FREE
# define KISS_FFT_FREE _mm_free
# endif
#else
# define KISS_FFT_ALIGN_CHECK(ptr)
# define KISS_FFT_ALIGN_SIZE_UP(size) (size)
# ifndef KISS_FFT_MALLOC
# define KISS_FFT_MALLOC malloc
# endif
# ifndef KISS_FFT_FREE
# define KISS_FFT_FREE free
# endif
#endif
#ifdef FIXED_POINT #ifdef FIXED_POINT
#include <sys/types.h> #include <stdint.h>
# if (FIXED_POINT == 32) # if (FIXED_POINT == 32)
# define kiss_fft_scalar int32_t # define kiss_fft_scalar int32_t
# else # else
@@ -79,7 +114,7 @@ typedef struct kiss_fft_state* kiss_fft_cfg;
* buffer size in *lenmem. * buffer size in *lenmem.
* */ * */
kiss_fft_cfg kiss_fft_alloc(int nfft,int inverse_fft,void * mem,size_t * lenmem); kiss_fft_cfg KISS_FFT_API kiss_fft_alloc(int nfft,int inverse_fft,void * mem,size_t * lenmem);
/* /*
* kiss_fft(cfg,in_out_buf) * kiss_fft(cfg,in_out_buf)
@@ -91,28 +126,32 @@ kiss_fft_cfg kiss_fft_alloc(int nfft,int inverse_fft,void * mem,size_t * lenmem)
* Note that each element is complex and can be accessed like * Note that each element is complex and can be accessed like
f[k].r and f[k].i f[k].r and f[k].i
* */ * */
void kiss_fft(kiss_fft_cfg cfg,const kiss_fft_cpx *fin,kiss_fft_cpx *fout); void KISS_FFT_API kiss_fft(kiss_fft_cfg cfg,const kiss_fft_cpx *fin,kiss_fft_cpx *fout);
/* /*
A more generic version of the above function. It reads its input from every Nth sample. A more generic version of the above function. It reads its input from every Nth sample.
* */ * */
void kiss_fft_stride(kiss_fft_cfg cfg,const kiss_fft_cpx *fin,kiss_fft_cpx *fout,int fin_stride); void KISS_FFT_API kiss_fft_stride(kiss_fft_cfg cfg,const kiss_fft_cpx *fin,kiss_fft_cpx *fout,int fin_stride);
/* If kiss_fft_alloc allocated a buffer, it is one contiguous /* If kiss_fft_alloc allocated a buffer, it is one contiguous
buffer and can be simply free()d when no longer needed*/ buffer and can be simply free()d when no longer needed*/
#define kiss_fft_free free #define kiss_fft_free KISS_FFT_FREE
/* /*
Cleans up some memory that gets managed internally. Not necessary to call, but it might clean up Cleans up some memory that gets managed internally. Not necessary to call, but it might clean up
your compiler output to call this before you exit. your compiler output to call this before you exit.
*/ */
void kiss_fft_cleanup(void); void KISS_FFT_API kiss_fft_cleanup(void);
/* /*
* Returns the smallest integer k, such that k>=n and k has only "fast" factors (2,3,5) * Returns the smallest integer k, such that k>=n and k has only "fast" factors (2,3,5)
*/ */
int kiss_fft_next_fast_size(int n); int KISS_FFT_API kiss_fft_next_fast_size(int n);
/* for real ffts, we need an even size */
#define kiss_fftr_next_fast_size_real(n) \
(kiss_fft_next_fast_size( ((n)+1)>>1)<<1)
#ifdef __cplusplus #ifdef __cplusplus
} }

View File

@@ -0,0 +1,36 @@
/*
* Copyright (c) 2003-2010, Mark Borgerding. All rights reserved.
* This file is part of KISS FFT - https://github.com/mborgerding/kissfft
*
* SPDX-License-Identifier: BSD-3-Clause
* See COPYING file for more information.
*/
#ifndef kiss_fft_log_h
#define kiss_fft_log_h
#define ERROR 1
#define WARNING 2
#define INFO 3
#define DEBUG 4
#define STRINGIFY(x) #x
#define TOSTRING(x) STRINGIFY(x)
#if defined(NDEBUG)
# define KISS_FFT_LOG_MSG(severity, ...) ((void)0)
#else
# define KISS_FFT_LOG_MSG(severity, ...) \
fprintf(stderr, "[" #severity "] " __FILE__ ":" TOSTRING(__LINE__) " "); \
fprintf(stderr, __VA_ARGS__); \
fprintf(stderr, "\n")
#endif
#define KISS_FFT_ERROR(...) KISS_FFT_LOG_MSG(ERROR, __VA_ARGS__)
#define KISS_FFT_WARNING(...) KISS_FFT_LOG_MSG(WARNING, __VA_ARGS__)
#define KISS_FFT_INFO(...) KISS_FFT_LOG_MSG(INFO, __VA_ARGS__)
#define KISS_FFT_DEBUG(...) KISS_FFT_LOG_MSG(DEBUG, __VA_ARGS__)
#endif /* kiss_fft_log_h */

View File

@@ -1,16 +1,10 @@
/* /*
Copyright (c) 2003-2004, Mark Borgerding * Copyright (c) 2003-2004, Mark Borgerding. All rights reserved.
* This file is part of KISS FFT - https://github.com/mborgerding/kissfft
All rights reserved. *
* SPDX-License-Identifier: BSD-3-Clause
Redistribution and use in source and binary forms, with or without modification, are permitted provided that the following conditions are met: * See COPYING file for more information.
*/
* Redistributions of source code must retain the above copyright notice, this list of conditions and the following disclaimer.
* Redistributions in binary form must reproduce the above copyright notice, this list of conditions and the following disclaimer in the documentation and/or other materials provided with the distribution.
* Neither the author nor the names of any contributors may be used to endorse or promote products derived from this software without specific prior written permission.
THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
*/
#include "kiss_fftr.h" #include "kiss_fftr.h"
#include "_kiss_fft_guts.h" #include "_kiss_fft_guts.h"
@@ -19,25 +13,27 @@ struct kiss_fftr_state{
kiss_fft_cfg substate; kiss_fft_cfg substate;
kiss_fft_cpx * tmpbuf; kiss_fft_cpx * tmpbuf;
kiss_fft_cpx * super_twiddles; kiss_fft_cpx * super_twiddles;
#ifdef USE_SIMD #ifdef USE_SIMD
long pad; void * pad;
#endif #endif
}; };
kiss_fftr_cfg kiss_fftr_alloc(int nfft,int inverse_fft,void * mem,size_t * lenmem) kiss_fftr_cfg kiss_fftr_alloc(int nfft,int inverse_fft,void * mem,size_t * lenmem)
{ {
KISS_FFT_ALIGN_CHECK(mem)
int i; int i;
kiss_fftr_cfg st = NULL; kiss_fftr_cfg st = NULL;
size_t subsize, memneeded; size_t subsize = 0, memneeded;
if (nfft & 1) { if (nfft & 1) {
fprintf(stderr,"Real FFT optimization must be even.\n"); KISS_FFT_ERROR("Real FFT optimization must be even.");
return NULL; return NULL;
} }
nfft >>= 1; nfft >>= 1;
kiss_fft_alloc (nfft, inverse_fft, NULL, &subsize); kiss_fft_alloc (nfft, inverse_fft, NULL, &subsize);
memneeded = sizeof(struct kiss_fftr_state) + subsize + sizeof(kiss_fft_cpx) * ( nfft * 2); memneeded = sizeof(struct kiss_fftr_state) + subsize + sizeof(kiss_fft_cpx) * ( nfft * 3 / 2);
if (lenmem == NULL) { if (lenmem == NULL) {
st = (kiss_fftr_cfg) KISS_FFT_MALLOC (memneeded); st = (kiss_fftr_cfg) KISS_FFT_MALLOC (memneeded);
@@ -54,9 +50,9 @@ kiss_fftr_cfg kiss_fftr_alloc(int nfft,int inverse_fft,void * mem,size_t * lenme
st->super_twiddles = st->tmpbuf + nfft; st->super_twiddles = st->tmpbuf + nfft;
kiss_fft_alloc(nfft, inverse_fft, st->substate, &subsize); kiss_fft_alloc(nfft, inverse_fft, st->substate, &subsize);
for (i = 0; i < nfft; ++i) { for (i = 0; i < nfft/2; ++i) {
double phase = double phase =
-3.14159265358979323846264338327 * ((double) i / nfft + .5); -3.14159265358979323846264338327 * ((double) (i+1) / nfft + .5);
if (inverse_fft) if (inverse_fft)
phase *= -1; phase *= -1;
kf_cexp (st->super_twiddles+i,phase); kf_cexp (st->super_twiddles+i,phase);
@@ -71,8 +67,8 @@ void kiss_fftr(kiss_fftr_cfg st,const kiss_fft_scalar *timedata,kiss_fft_cpx *fr
kiss_fft_cpx fpnk,fpk,f1k,f2k,tw,tdc; kiss_fft_cpx fpnk,fpk,f1k,f2k,tw,tdc;
if ( st->substate->inverse) { if ( st->substate->inverse) {
fprintf(stderr,"kiss fft usage error: improper alloc\n"); KISS_FFT_ERROR("kiss fft usage error: improper alloc");
exit(1); return;/* The caller did not call the correct function */
} }
ncfft = st->substate->nfft; ncfft = st->substate->nfft;
@@ -83,12 +79,12 @@ void kiss_fftr(kiss_fftr_cfg st,const kiss_fft_scalar *timedata,kiss_fft_cpx *fr
* contains the sum of the even-numbered elements of the input time sequence * contains the sum of the even-numbered elements of the input time sequence
* The imag part is the sum of the odd-numbered elements * The imag part is the sum of the odd-numbered elements
* *
* The sum of tdc.r and tdc.i is the sum of the input time sequence. * The sum of tdc.r and tdc.i is the sum of the input time sequence.
* yielding DC of input time sequence * yielding DC of input time sequence
* The difference of tdc.r - tdc.i is the sum of the input (dot product) [1,-1,1,-1... * The difference of tdc.r - tdc.i is the sum of the input (dot product) [1,-1,1,-1...
* yielding Nyquist bin of input time sequence * yielding Nyquist bin of input time sequence
*/ */
tdc.r = st->tmpbuf[0].r; tdc.r = st->tmpbuf[0].r;
tdc.i = st->tmpbuf[0].i; tdc.i = st->tmpbuf[0].i;
C_FIXDIV(tdc,2); C_FIXDIV(tdc,2);
@@ -96,14 +92,14 @@ void kiss_fftr(kiss_fftr_cfg st,const kiss_fft_scalar *timedata,kiss_fft_cpx *fr
CHECK_OVERFLOW_OP(tdc.r ,-, tdc.i); CHECK_OVERFLOW_OP(tdc.r ,-, tdc.i);
freqdata[0].r = tdc.r + tdc.i; freqdata[0].r = tdc.r + tdc.i;
freqdata[ncfft].r = tdc.r - tdc.i; freqdata[ncfft].r = tdc.r - tdc.i;
#ifdef USE_SIMD #ifdef USE_SIMD
freqdata[ncfft].i = freqdata[0].i = _mm_set1_ps(0); freqdata[ncfft].i = freqdata[0].i = _mm_set1_ps(0);
#else #else
freqdata[ncfft].i = freqdata[0].i = 0; freqdata[ncfft].i = freqdata[0].i = 0;
#endif #endif
for ( k=1;k <= ncfft/2 ; ++k ) { for ( k=1;k <= ncfft/2 ; ++k ) {
fpk = st->tmpbuf[k]; fpk = st->tmpbuf[k];
fpnk.r = st->tmpbuf[ncfft-k].r; fpnk.r = st->tmpbuf[ncfft-k].r;
fpnk.i = - st->tmpbuf[ncfft-k].i; fpnk.i = - st->tmpbuf[ncfft-k].i;
C_FIXDIV(fpk,2); C_FIXDIV(fpk,2);
@@ -111,7 +107,7 @@ void kiss_fftr(kiss_fftr_cfg st,const kiss_fft_scalar *timedata,kiss_fft_cpx *fr
C_ADD( f1k, fpk , fpnk ); C_ADD( f1k, fpk , fpnk );
C_SUB( f2k, fpk , fpnk ); C_SUB( f2k, fpk , fpnk );
C_MUL( tw , f2k , st->super_twiddles[k]); C_MUL( tw , f2k , st->super_twiddles[k-1]);
freqdata[k].r = HALF_OF(f1k.r + tw.r); freqdata[k].r = HALF_OF(f1k.r + tw.r);
freqdata[k].i = HALF_OF(f1k.i + tw.i); freqdata[k].i = HALF_OF(f1k.i + tw.i);
@@ -126,8 +122,8 @@ void kiss_fftri(kiss_fftr_cfg st,const kiss_fft_cpx *freqdata,kiss_fft_scalar *t
int k, ncfft; int k, ncfft;
if (st->substate->inverse == 0) { if (st->substate->inverse == 0) {
fprintf (stderr, "kiss fft usage error: improper alloc\n"); KISS_FFT_ERROR("kiss fft usage error: improper alloc");
exit (1); return;/* The caller did not call the correct function */
} }
ncfft = st->substate->nfft; ncfft = st->substate->nfft;
@@ -146,10 +142,10 @@ void kiss_fftri(kiss_fftr_cfg st,const kiss_fft_cpx *freqdata,kiss_fft_scalar *t
C_ADD (fek, fk, fnkc); C_ADD (fek, fk, fnkc);
C_SUB (tmp, fk, fnkc); C_SUB (tmp, fk, fnkc);
C_MUL (fok, tmp, st->super_twiddles[k]); C_MUL (fok, tmp, st->super_twiddles[k-1]);
C_ADD (st->tmpbuf[k], fek, fok); C_ADD (st->tmpbuf[k], fek, fok);
C_SUB (st->tmpbuf[ncfft - k], fek, fok); C_SUB (st->tmpbuf[ncfft - k], fek, fok);
#ifdef USE_SIMD #ifdef USE_SIMD
st->tmpbuf[ncfft - k].i *= _mm_set1_ps(-1.0); st->tmpbuf[ncfft - k].i *= _mm_set1_ps(-1.0);
#else #else
st->tmpbuf[ncfft - k].i *= -1; st->tmpbuf[ncfft - k].i *= -1;

View File

@@ -1,3 +1,11 @@
/*
* Copyright (c) 2003-2004, Mark Borgerding. All rights reserved.
* This file is part of KISS FFT - https://github.com/mborgerding/kissfft
*
* SPDX-License-Identifier: BSD-3-Clause
* See COPYING file for more information.
*/
#ifndef KISS_FTR_H #ifndef KISS_FTR_H
#define KISS_FTR_H #define KISS_FTR_H
@@ -18,7 +26,7 @@ extern "C" {
typedef struct kiss_fftr_state *kiss_fftr_cfg; typedef struct kiss_fftr_state *kiss_fftr_cfg;
kiss_fftr_cfg kiss_fftr_alloc(int nfft,int inverse_fft,void * mem, size_t * lenmem); kiss_fftr_cfg KISS_FFT_API kiss_fftr_alloc(int nfft,int inverse_fft,void * mem, size_t * lenmem);
/* /*
nfft must be even nfft must be even
@@ -26,19 +34,19 @@ kiss_fftr_cfg kiss_fftr_alloc(int nfft,int inverse_fft,void * mem, size_t * lenm
*/ */
void kiss_fftr(kiss_fftr_cfg cfg,const kiss_fft_scalar *timedata,kiss_fft_cpx *freqdata); void KISS_FFT_API kiss_fftr(kiss_fftr_cfg cfg,const kiss_fft_scalar *timedata,kiss_fft_cpx *freqdata);
/* /*
input timedata has nfft scalar points input timedata has nfft scalar points
output freqdata has nfft/2+1 complex points output freqdata has nfft/2+1 complex points
*/ */
void kiss_fftri(kiss_fftr_cfg cfg,const kiss_fft_cpx *freqdata,kiss_fft_scalar *timedata); void KISS_FFT_API kiss_fftri(kiss_fftr_cfg cfg,const kiss_fft_cpx *freqdata,kiss_fft_scalar *timedata);
/* /*
input freqdata has nfft/2+1 complex points input freqdata has nfft/2+1 complex points
output timedata has nfft scalar points output timedata has nfft scalar points
*/ */
#define kiss_fftr_free free #define kiss_fftr_free KISS_FFT_FREE
#ifdef __cplusplus #ifdef __cplusplus
} }

View File

@@ -246,6 +246,63 @@ void v_cartesian_to_polar_interleaved_inplace(T *const R__ srcdst,
} }
} }
template<typename S, typename T> // S source, T target
void v_cartesian_to_magnitudes(T *const R__ mag,
const S *const R__ real,
const S *const R__ imag,
const int count)
{
for (int i = 0; i < count; ++i) {
mag[i] = T(sqrt(real[i] * real[i] + imag[i] * imag[i]));
}
}
template<typename S, typename T> // S source, T target
void v_cartesian_interleaved_to_magnitudes(T *const R__ mag,
const S *const R__ src,
const int count)
{
for (int i = 0; i < count; ++i) {
mag[i] = T(sqrt(src[i*2] * src[i*2] + src[i*2+1] * src[i*2+1]));
}
}
#ifdef HAVE_IPP
template<>
inline void v_cartesian_to_magnitudes(float *const R__ mag,
const float *const R__ real,
const float *const R__ imag,
const int count)
{
ippsMagnitude_32f(real, imag, mag, count);
}
template<>
inline void v_cartesian_to_magnitudes(double *const R__ mag,
const double *const R__ real,
const double *const R__ imag,
const int count)
{
ippsMagnitude_64f(real, imag, mag, count);
}
template<>
inline void v_cartesian_interleaved_to_magnitudes(float *const R__ mag,
const float *const R__ src,
const int count)
{
ippsMagnitude_32fc((const Ipp32fc *)src, mag, count);
}
template<>
inline void v_cartesian_interleaved_to_magnitudes(double *const R__ mag,
const double *const R__ src,
const int count)
{
ippsMagnitude_64fc((const Ipp64fc *)src, mag, count);
}
#endif
} }
#endif #endif