jkl1337/ecgsyn.js
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

add pffft

Browse Source
This commit is contained in:
John K. Luebs
2024-11-09 14:57:18 -06:00
parent 78a00f71cc
commit a1790b8977
69 changed files with 25719 additions and 0 deletions
Download Patch File Download Diff File Expand all files Collapse all files

377
pffft/test_pffft.cpp Normal file
View File

@@ -0,0 +1,377 @@
/*
Copyright (c) 2013 Julien Pommier ( pommier@modartt.com )
Copyright (c) 2020 Dario Mambro ( dario.mambro@gmail.com )
Copyright (c) 2020 Hayati Ayguen ( h_ayguen@web.de )
Small test & bench for PFFFT, comparing its performance with the scalar
FFTPACK, FFTW, and Apple vDSP
How to build:
on linux, with fftw3:
gcc -o test_pffft -DHAVE_FFTW -msse -mfpmath=sse -O3 -Wall -W pffft.c
test_pffft.c fftpack.c -L/usr/local/lib -I/usr/local/include/ -lfftw3f -lm
on macos, without fftw3:
clang -o test_pffft -DHAVE_VECLIB -O3 -Wall -W pffft.c test_pffft.c fftpack.c
-L/usr/local/lib -I/usr/local/include/ -framework Accelerate
on macos, with fftw3:
clang -o test_pffft -DHAVE_FFTW -DHAVE_VECLIB -O3 -Wall -W pffft.c
test_pffft.c fftpack.c -L/usr/local/lib -I/usr/local/include/ -lfftw3f
-framework Accelerate
as alternative: replace clang by gcc.
on windows, with visual c++:
cl /Ox -D_USE_MATH_DEFINES /arch:SSE test_pffft.c pffft.c fftpack.c
build without SIMD instructions:
gcc -o test_pffft -DPFFFT_SIMD_DISABLE -O3 -Wall -W pffft.c test_pffft.c
fftpack.c -lm
*/
#include "pffft.hpp"
#include <assert.h>
#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <time.h>
/* define own constants required to turn off g++ extensions .. */
#ifndef M_PI
#define M_PI 3.14159265358979323846 /* pi */
#endif
/* maximum allowed phase error in degree */
#define DEG_ERR_LIMIT 1E-4
/* maximum allowed magnitude error in amplitude (of 1.0 or 1.1) */
#define MAG_ERR_LIMIT 1E-6
#define PRINT_SPEC 0
#define PWR2LOG(PWR) ((PWR) < 1E-30 ? 10.0 * log10(1E-30) : 10.0 * log10(PWR))
template<typename T>
bool
Ttest(int N, bool useOrdered)
{
typedef pffft::Fft<T> Fft;
typedef typename pffft::Fft<T>::Scalar FftScalar;
typedef typename Fft::Complex FftComplex;
const bool cplx = pffft::Fft<T>::isComplexTransform();
const double EXPECTED_DYN_RANGE = Fft::isDoubleScalar() ? 215.0 : 140.0;
assert(Fft::isPowerOfTwo(N));
Fft fft = Fft(N); // instantiate and prepareLength() for length N
#if __cplusplus >= 201103L || (defined(_MSC_VER) && _MSC_VER >= 1900)
// possible ways to declare/instatiate aligned vectors with C++11
// some lines require a typedef of above
auto X = fft.valueVector(); // for X = input vector
pffft::AlignedVector<typename Fft::Complex> Y = fft.spectrumVector(); // for Y = forward(X)
pffft::AlignedVector<FftScalar> R = fft.internalLayoutVector(); // for R = forwardInternalLayout(X)
pffft::AlignedVector<T> Z = fft.valueVector(); // for Z = inverse(Y) = inverse( forward(X) )
// or Z = inverseInternalLayout(R)
#else
// possible ways to declare/instatiate aligned vectors with C++98
pffft::AlignedVector<T> X = fft.valueVector(); // for X = input vector
pffft::AlignedVector<FftComplex> Y = fft.spectrumVector(); // for Y = forward(X)
pffft::AlignedVector<typename Fft::Scalar> R = fft.internalLayoutVector(); // for R = forwardInternalLayout(X)
pffft::AlignedVector<T> Z = fft.valueVector(); // for Z = inverse(Y) = inverse( forward(X) )
// or Z = inverseInternalLayout(R)
#endif
// work with complex - without the capabilities of a higher c++ standard
FftScalar* Xs = reinterpret_cast<FftScalar*>(X.data()); // for X = input vector
FftScalar* Ys = reinterpret_cast<FftScalar*>(Y.data()); // for Y = forward(X)
FftScalar* Zs = reinterpret_cast<FftScalar*>(Z.data()); // for Z = inverse(Y) = inverse( forward(X) )
int k, j, m, iter, kmaxOther;
bool retError = false;
double freq, dPhi, phi, phi0;
double pwr, pwrCar, pwrOther, err, errSum, mag, expextedMag;
double amp = 1.0;
for (k = m = 0; k < (cplx ? N : (1 + N / 2)); k += N / 16, ++m) {
amp = ((m % 3) == 0) ? 1.0F : 1.1F;
freq = (k < N / 2) ? ((double)k / N) : ((double)(k - N) / N);
dPhi = 2.0 * M_PI * freq;
if (dPhi < 0.0)
dPhi += 2.0 * M_PI;
iter = -1;
while (1) {
++iter;
if (iter)
printf("bin %d: dphi = %f for freq %f\n", k, dPhi, freq);
/* generate cosine carrier as time signal - start at defined phase phi0 */
phi = phi0 =
(m % 4) * 0.125 * M_PI; /* have phi0 < 90 deg to be normalized */
for (j = 0; j < N; ++j) {
if (cplx) {
Xs[2 * j] = (FftScalar)( amp * cos(phi) ); /* real part */
Xs[2 * j + 1] = (FftScalar)( amp * sin(phi) ); /* imag part */
} else
Xs[j] = (FftScalar)( amp * cos(phi) ); /* only real part */
/* phase increment .. stay normalized - cos()/sin() might degrade! */
phi += dPhi;
if (phi >= M_PI)
phi -= 2.0 * M_PI;
}
/* forward transform from X --> Y .. using work buffer W */
if (useOrdered)
fft.forward(X, Y);
else {
fft.forwardToInternalLayout(X, R); /* use R for reordering */
fft.reorderSpectrum(R, Y); /* have canonical order in Y[] for power calculations */
}
pwrOther = -1.0;
pwrCar = 0;
/* for positive frequencies: 0 to 0.5 * samplerate */
/* and also for negative frequencies: -0.5 * samplerate to 0 */
for (j = 0; j < (cplx ? N : (1 + N / 2)); ++j) {
if (!cplx && !j) /* special treatment for DC for real input */
pwr = Ys[j] * Ys[j];
else if (!cplx && j == N / 2) /* treat 0.5 * samplerate */
pwr = Ys[1] *
Ys[1]; /* despite j (for freq calculation) we have index 1 */
else
pwr = Ys[2 * j] * Ys[2 * j] + Ys[2 * j + 1] * Ys[2 * j + 1];
if (iter || PRINT_SPEC)
printf("%s fft %d: pwr[j = %d] = %g == %f dB\n",
(cplx ? "cplx" : "real"),
N,
j,
pwr,
PWR2LOG(pwr));
if (k == j)
pwrCar = pwr;
else if (pwr > pwrOther) {
pwrOther = pwr;
kmaxOther = j;
}
}
if (PWR2LOG(pwrCar) - PWR2LOG(pwrOther) < EXPECTED_DYN_RANGE) {
printf("%s fft %d amp %f iter %d:\n",
(cplx ? "cplx" : "real"),
N,
amp,
iter);
printf(" carrier power at bin %d: %g == %f dB\n",
k,
pwrCar,
PWR2LOG(pwrCar));
printf(" carrier mag || at bin %d: %g\n", k, sqrt(pwrCar));
printf(" max other pwr at bin %d: %g == %f dB\n",
kmaxOther,
pwrOther,
PWR2LOG(pwrOther));
printf(" dynamic range: %f dB\n\n",
PWR2LOG(pwrCar) - PWR2LOG(pwrOther));
retError = true;
if (iter == 0)
continue;
}
if (k > 0 && k != N / 2) {
phi = atan2(Ys[2 * k + 1], Ys[2 * k]);
if (fabs(phi - phi0) > DEG_ERR_LIMIT * M_PI / 180.0) {
retError = true;
printf("%s fft %d bin %d amp %f : phase mismatch! phase = %f deg "
"expected = %f deg\n",
(cplx ? "cplx" : "real"),
N,
k,
amp,
phi * 180.0 / M_PI,
phi0 * 180.0 / M_PI);
}
}
expextedMag = cplx ? amp : ((k == 0 || k == N / 2) ? amp : (amp / 2));
mag = sqrt(pwrCar) / N;
if (fabs(mag - expextedMag) > MAG_ERR_LIMIT) {
retError = true;
printf("%s fft %d bin %d amp %f : mag = %g expected = %g\n",
(cplx ? "cplx" : "real"),
N,
k,
amp,
mag,
expextedMag);
}
/* now convert spectrum back */
if (useOrdered)
fft.inverse(Y, Z);
else
fft.inverseFromInternalLayout(R, Z); /* inverse() from internal Layout */
errSum = 0.0;
for (j = 0; j < (cplx ? (2 * N) : N); ++j) {
/* scale back */
Zs[j] /= N;
/* square sum errors over real (and imag parts) */
err = (Xs[j] - Zs[j]) * (Xs[j] - Zs[j]);
errSum += err;
}
if (errSum > N * 1E-7) {
retError = true;
printf("%s fft %d bin %d : inverse FFT doesn't match original signal! "
"errSum = %g ; mean err = %g\n",
(cplx ? "cplx" : "real"),
N,
k,
errSum,
errSum / N);
}
break;
}
}
// using the std::vector<> base classes .. no need for alignedFree() for X, Y, Z and R
return retError;
}
bool
test(int N, bool useComplex, bool useOrdered)
{
if (useComplex) {
return
#ifdef PFFFT_ENABLE_FLOAT
Ttest< std::complex<float> >(N, useOrdered)
#endif
#if defined(PFFFT_ENABLE_FLOAT) && defined(PFFFT_ENABLE_DOUBLE)
&&
#endif
#ifdef PFFFT_ENABLE_DOUBLE
Ttest< std::complex<double> >(N, useOrdered)
#endif
;
} else {
return
#ifdef PFFFT_ENABLE_FLOAT
Ttest<float>(N, useOrdered)
#endif
#if defined(PFFFT_ENABLE_FLOAT) && defined(PFFFT_ENABLE_DOUBLE)
&&
#endif
#ifdef PFFFT_ENABLE_DOUBLE
Ttest<double>(N, useOrdered)
#endif
;
}
}
int
main(int argc, char** argv)
{
int N, result, resN, resAll, k, resNextPw2, resIsPw2, resFFT;
int inp_power_of_two[] = { 1, 2, 3, 4, 5, 6, 7, 8, 9, 511, 512, 513 };
int ref_power_of_two[] = { 1, 2, 4, 4, 8, 8, 8, 8, 16, 512, 512, 1024 };
resNextPw2 = 0;
resIsPw2 = 0;
for (k = 0; k < (sizeof(inp_power_of_two) / sizeof(inp_power_of_two[0]));
++k) {
#ifdef PFFFT_ENABLE_FLOAT
N = pffft::Fft<float>::nextPowerOfTwo(inp_power_of_two[k]);
#else
N = pffft::Fft<double>::nextPowerOfTwo(inp_power_of_two[k]);
#endif
if (N != ref_power_of_two[k]) {
resNextPw2 = 1;
printf("pffft_next_power_of_two(%d) does deliver %d, which is not "
"reference result %d!\n",
inp_power_of_two[k],
N,
ref_power_of_two[k]);
}
#ifdef PFFFT_ENABLE_FLOAT
result = pffft::Fft<float>::isPowerOfTwo(inp_power_of_two[k]);
#else
result = pffft::Fft<double>::isPowerOfTwo(inp_power_of_two[k]);
#endif
if (inp_power_of_two[k] == ref_power_of_two[k]) {
if (!result) {
resIsPw2 = 1;
printf("pffft_is_power_of_two(%d) delivers false; expected true!\n",
inp_power_of_two[k]);
}
} else {
if (result) {
resIsPw2 = 1;
printf("pffft_is_power_of_two(%d) delivers true; expected false!\n",
inp_power_of_two[k]);
}
}
}
if (!resNextPw2)
printf("tests for pffft_next_power_of_two() succeeded successfully.\n");
if (!resIsPw2)
printf("tests for pffft_is_power_of_two() succeeded successfully.\n");
resFFT = 0;
for (N = 32; N <= 65536; N *= 2) {
result = test(N, 1 /* cplx fft */, 1 /* useOrdered */);
resN = result;
resFFT |= result;
result = test(N, 0 /* cplx fft */, 1 /* useOrdered */);
resN |= result;
resFFT |= result;
result = test(N, 1 /* cplx fft */, 0 /* useOrdered */);
resN |= result;
resFFT |= result;
result = test(N, 0 /* cplx fft */, 0 /* useOrdered */);
resN |= result;
resFFT |= result;
if (!resN)
printf("tests for size %d succeeded successfully.\n", N);
}
if (!resFFT)
printf("all pffft transform tests (FORWARD/BACKWARD, REAL/COMPLEX, "
#ifdef PFFFT_ENABLE_FLOAT
"float"
#endif
#if defined(PFFFT_ENABLE_FLOAT) && defined(PFFFT_ENABLE_DOUBLE)
"/"
#endif
#ifdef PFFFT_ENABLE_DOUBLE
"double"
#endif
") succeeded successfully.\n");
resAll = resNextPw2 | resIsPw2 | resFFT;
if (!resAll)
printf("all tests succeeded successfully.\n");
else
printf("there are failed tests!\n");
return resAll;
}