Dima353/openal-soft
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

Use standard complex types instead of custom

Browse Source
This commit is contained in:
Chris Robinson
2018-11-17 04:14:57 -08:00
parent 09943683b5
commit ccdaca80c9
6 changed files with 145 additions and 217 deletions
Download Patch File Download Diff File Expand all files Collapse all files
Alc/effects/fshifter.cpp
+21 -20
View File
@@ -22,6 +22,8 @@
#include <cmath>
#include <cstdlib>
#include <complex>
#include <algorithm>
#include "alMain.h"
#include "alAuxEffectSlot.h"
@@ -33,6 +35,8 @@
namespace {
using complex_d = std::complex<double>;
#define HIL_SIZE 1024
#define OVERSAMP (1<<2)
@@ -64,10 +68,10 @@ struct ALfshifterState final : public ALeffectState {
/*Effects buffers*/
ALfloat InFIFO[HIL_SIZE];
ALcomplex OutFIFO[HIL_SIZE];
ALcomplex OutputAccum[HIL_SIZE];
ALcomplex Analytic[HIL_SIZE];
ALcomplex Outdata[BUFFERSIZE];
complex_d OutFIFO[HIL_SIZE];
complex_d OutputAccum[HIL_SIZE];
complex_d Analytic[HIL_SIZE];
complex_d Outdata[BUFFERSIZE];
alignas(16) ALfloat BufferOut[BUFFERSIZE];
@@ -105,13 +109,13 @@ ALboolean ALfshifterState_deviceUpdate(ALfshifterState *state, ALCdevice *UNUSED
state->Phase = 0;
state->ld_sign = 1.0;
memset(state->InFIFO, 0, sizeof(state->InFIFO));
memset(state->OutFIFO, 0, sizeof(state->OutFIFO));
memset(state->OutputAccum, 0, sizeof(state->OutputAccum));
memset(state->Analytic, 0, sizeof(state->Analytic));
std::fill(std::begin(state->InFIFO), std::end(state->InFIFO), 0.0f);
std::fill(std::begin(state->OutFIFO), std::end(state->OutFIFO), complex_d{});
std::fill(std::begin(state->OutputAccum), std::end(state->OutputAccum), complex_d{});
std::fill(std::begin(state->Analytic), std::end(state->Analytic), complex_d{});
memset(state->CurrentGains, 0, sizeof(state->CurrentGains));
memset(state->TargetGains, 0, sizeof(state->TargetGains));
std::fill(std::begin(state->CurrentGains), std::end(state->CurrentGains), 0.0f);
std::fill(std::begin(state->TargetGains), std::end(state->TargetGains), 0.0f);
return AL_TRUE;
}
@@ -147,7 +151,7 @@ ALvoid ALfshifterState_update(ALfshifterState *state, const ALCcontext *context,
ALvoid ALfshifterState_process(ALfshifterState *state, ALsizei SamplesToDo, const ALfloat (*RESTRICT SamplesIn)[BUFFERSIZE], ALfloat (*RESTRICT SamplesOut)[BUFFERSIZE], ALsizei NumChannels)
{
static const ALcomplex complex_zero = { 0.0, 0.0 };
static const complex_d complex_zero{0.0, 0.0};
ALfloat *RESTRICT BufferOut = state->BufferOut;
ALsizei j, k, base;
@@ -175,8 +179,8 @@ ALvoid ALfshifterState_process(ALfshifterState *state, ALsizei SamplesToDo, cons
/* Real signal windowing and store in Analytic buffer */
for(k = 0;k < HIL_SIZE;k++)
{
state->Analytic[k].Real = state->InFIFO[k] * HannWindow[k];
state->Analytic[k].Imag = 0.0;
state->Analytic[k].real(state->InFIFO[k] * HannWindow[k]);
state->Analytic[k].imag(0.0);
}
/* Processing signal by Discrete Hilbert Transform (analytical signal). */
@@ -184,10 +188,7 @@ ALvoid ALfshifterState_process(ALfshifterState *state, ALsizei SamplesToDo, cons
/* Windowing and add to output accumulator */
for(k = 0;k < HIL_SIZE;k++)
{
state->OutputAccum[k].Real += 2.0/OVERSAMP*HannWindow[k]*state->Analytic[k].Real;
state->OutputAccum[k].Imag += 2.0/OVERSAMP*HannWindow[k]*state->Analytic[k].Imag;
}
state->OutputAccum[k] += 2.0/OVERSAMP*HannWindow[k]*state->Analytic[k];
/* Shift accumulator, input & output FIFO */
for(k = 0;k < HIL_STEP;k++) state->OutFIFO[k] = state->OutputAccum[k];
@@ -200,9 +201,9 @@ ALvoid ALfshifterState_process(ALfshifterState *state, ALsizei SamplesToDo, cons
/* Process frequency shifter using the analytic signal obtained. */
for(k = 0;k < SamplesToDo;k++)
{
ALdouble phase = state->Phase * ((1.0/FRACTIONONE) * 2.0*M_PI);
BufferOut[k] = (ALfloat)(state->Outdata[k].Real*cos(phase) +
state->Outdata[k].Imag*sin(phase)*state->ld_sign);
double phase = state->Phase * ((1.0/FRACTIONONE) * 2.0*M_PI);
BufferOut[k] = (float)(state->Outdata[k].real()*std::cos(phase) +
state->Outdata[k].imag()*std::sin(phase)*state->ld_sign);
state->Phase += state->PhaseStep;
state->Phase &= FRACTIONMASK;
Alc/effects/pshifter.cpp
+44 -57
View File
@@ -22,6 +22,8 @@
#include <cmath>
#include <cstdlib>
#include <complex>
#include <algorithm>
#include "alMain.h"
#include "alAuxEffectSlot.h"
@@ -34,6 +36,8 @@
namespace {
using complex_d = std::complex<double>;
#define STFT_SIZE 1024
#define STFT_HALF_SIZE (STFT_SIZE>>1)
#define OVERSAMP (1<<2)
@@ -41,7 +45,7 @@ namespace {
#define STFT_STEP (STFT_SIZE / OVERSAMP)
#define FIFO_LATENCY (STFT_STEP * (OVERSAMP-1))
inline ALint double2int(ALdouble d)
inline int double2int(double d)
{
#if ((defined(__GNUC__) || defined(__clang__)) && (defined(__i386__) || defined(__x86_64__)) && \
!defined(__SSE2_MATH__)) || (defined(_MSC_VER) && defined(_M_IX86_FP) && _M_IX86_FP < 2)
@@ -98,28 +102,18 @@ struct ALfrequencyDomain {
};
/* Converts ALcomplex to ALphasor */
inline ALphasor rect2polar(ALcomplex number)
/* Converts complex to ALphasor */
inline ALphasor rect2polar(const complex_d &number)
{
ALphasor polar;
polar.Amplitude = std::sqrt(number.Real*number.Real + number.Imag*number.Imag);
polar.Phase = std::atan2(number.Imag, number.Real);
polar.Amplitude = std::abs(number);
polar.Phase = std::arg(number);
return polar;
}
/* Converts ALphasor to ALcomplex */
inline ALcomplex polar2rect(ALphasor number)
{
ALcomplex cartesian;
cartesian.Real = number.Amplitude * std::cos(number.Phase);
cartesian.Imag = number.Amplitude * std::sin(number.Phase);
return cartesian;
}
/* Converts ALphasor to complex */
inline complex_d polar2rect(const ALphasor &number)
{ return std::polar<double>(number.Amplitude, number.Phase); }
struct ALpshifterState final : public ALeffectState {
@@ -136,7 +130,7 @@ struct ALpshifterState final : public ALeffectState {
ALdouble SumPhase[STFT_HALF_SIZE+1];
ALdouble OutputAccum[STFT_SIZE];
ALcomplex FFTbuffer[STFT_SIZE];
complex_d FFTbuffer[STFT_SIZE];
ALfrequencyDomain Analysis_buffer[STFT_HALF_SIZE+1];
ALfrequencyDomain Syntesis_buffer[STFT_HALF_SIZE+1];
@@ -177,17 +171,17 @@ ALboolean ALpshifterState_deviceUpdate(ALpshifterState *state, ALCdevice *device
state->PitchShift = 1.0f;
state->FreqPerBin = device->Frequency / (ALfloat)STFT_SIZE;
memset(state->InFIFO, 0, sizeof(state->InFIFO));
memset(state->OutFIFO, 0, sizeof(state->OutFIFO));
memset(state->FFTbuffer, 0, sizeof(state->FFTbuffer));
memset(state->LastPhase, 0, sizeof(state->LastPhase));
memset(state->SumPhase, 0, sizeof(state->SumPhase));
memset(state->OutputAccum, 0, sizeof(state->OutputAccum));
memset(state->Analysis_buffer, 0, sizeof(state->Analysis_buffer));
memset(state->Syntesis_buffer, 0, sizeof(state->Syntesis_buffer));
std::fill(std::begin(state->InFIFO), std::end(state->InFIFO), 0.0f);
std::fill(std::begin(state->OutFIFO), std::end(state->OutFIFO), 0.0f);
std::fill(std::begin(state->LastPhase), std::end(state->LastPhase), 0.0);
std::fill(std::begin(state->SumPhase), std::end(state->SumPhase), 0.0);
std::fill(std::begin(state->OutputAccum), std::end(state->OutputAccum), 0.0);
std::fill(std::begin(state->FFTbuffer), std::end(state->FFTbuffer), complex_d{});
std::fill(std::begin(state->Analysis_buffer), std::end(state->Analysis_buffer), ALfrequencyDomain{});
std::fill(std::begin(state->Syntesis_buffer), std::end(state->Syntesis_buffer), ALfrequencyDomain{});
memset(state->CurrentGains, 0, sizeof(state->CurrentGains));
memset(state->TargetGains, 0, sizeof(state->TargetGains));
std::fill(std::begin(state->CurrentGains), std::end(state->CurrentGains), 0.0f);
std::fill(std::begin(state->TargetGains), std::end(state->TargetGains), 0.0f);
return AL_TRUE;
}
@@ -214,13 +208,12 @@ ALvoid ALpshifterState_process(ALpshifterState *state, ALsizei SamplesToDo, cons
* http://blogs.zynaptiq.com/bernsee/pitch-shifting-using-the-ft/
*/
static const ALdouble expected = M_PI*2.0 / OVERSAMP;
const ALdouble freq_per_bin = state->FreqPerBin;
ALfloat *RESTRICT bufferOut = state->BufferOut;
ALsizei count = state->count;
ALsizei i, j, k;
static constexpr ALdouble expected{M_PI*2.0 / OVERSAMP};
const ALdouble freq_per_bin{state->FreqPerBin};
ALfloat *RESTRICT bufferOut{state->BufferOut};
ALsizei count{state->count};
for(i = 0;i < SamplesToDo;)
for(ALsizei i{0};i < SamplesToDo;)
{
do {
/* Fill FIFO buffer with samples data */
@@ -235,10 +228,10 @@ ALvoid ALpshifterState_process(ALpshifterState *state, ALsizei SamplesToDo, cons
count = FIFO_LATENCY;
/* Real signal windowing and store in FFTbuffer */
for(k = 0;k < STFT_SIZE;k++)
for(ALsizei k{0};k < STFT_SIZE;k++)
{
state->FFTbuffer[k].Real = state->InFIFO[k] * HannWindow[k];
state->FFTbuffer[k].Imag = 0.0;
state->FFTbuffer[k].real(state->InFIFO[k] * HannWindow[k]);
state->FFTbuffer[k].imag(0.0);
}
/* ANALYSIS */
@@ -248,20 +241,16 @@ ALvoid ALpshifterState_process(ALpshifterState *state, ALsizei SamplesToDo, cons
/* Analyze the obtained data. Since the real FFT is symmetric, only
* STFT_HALF_SIZE+1 samples are needed.
*/
for(k = 0;k < STFT_HALF_SIZE+1;k++)
for(ALsizei k{0};k < STFT_HALF_SIZE+1;k++)
{
ALphasor component;
ALdouble tmp;
ALint qpd;
/* Compute amplitude and phase */
component = rect2polar(state->FFTbuffer[k]);
ALphasor component{rect2polar(state->FFTbuffer[k])};
/* Compute phase difference and subtract expected phase difference */
tmp = (component.Phase - state->LastPhase[k]) - k*expected;
double tmp{(component.Phase - state->LastPhase[k]) - k*expected};
/* Map delta phase into +/- Pi interval */
qpd = double2int(tmp / M_PI);
int qpd{double2int(tmp / M_PI)};
tmp -= M_PI * (qpd + (qpd%2));
/* Get deviation from bin frequency from the +/- Pi interval */
@@ -280,15 +269,15 @@ ALvoid ALpshifterState_process(ALpshifterState *state, ALsizei SamplesToDo, cons
/* PROCESSING */
/* pitch shifting */
for(k = 0;k < STFT_HALF_SIZE+1;k++)
for(ALsizei k{0};k < STFT_HALF_SIZE+1;k++)
{
state->Syntesis_buffer[k].Amplitude = 0.0;
state->Syntesis_buffer[k].Frequency = 0.0;
}
for(k = 0;k < STFT_HALF_SIZE+1;k++)
for(ALsizei k{0};k < STFT_HALF_SIZE+1;k++)
{
j = (k*state->PitchShiftI) >> FRACTIONBITS;
ALsizei j{(k*state->PitchShiftI) >> FRACTIONBITS};
if(j >= STFT_HALF_SIZE+1) break;
state->Syntesis_buffer[j].Amplitude += state->Analysis_buffer[k].Amplitude;
@@ -298,7 +287,7 @@ ALvoid ALpshifterState_process(ALpshifterState *state, ALsizei SamplesToDo, cons
/* SYNTHESIS */
/* Synthesis the processing data */
for(k = 0;k < STFT_HALF_SIZE+1;k++)
for(ALsizei k{0};k < STFT_HALF_SIZE+1;k++)
{
ALphasor component;
ALdouble tmp;
@@ -316,21 +305,19 @@ ALvoid ALpshifterState_process(ALpshifterState *state, ALsizei SamplesToDo, cons
state->FFTbuffer[k] = polar2rect(component);
}
/* zero negative frequencies for recontruct a real signal */
for(k = STFT_HALF_SIZE+1;k < STFT_SIZE;k++)
{
state->FFTbuffer[k].Real = 0.0;
state->FFTbuffer[k].Imag = 0.0;
}
for(ALsizei k{STFT_HALF_SIZE+1};k < STFT_SIZE;k++)
state->FFTbuffer[k] = complex_d{};
/* Apply iFFT to buffer data */
complex_fft(state->FFTbuffer, STFT_SIZE, 1.0);
/* Windowing and add to output */
for(k = 0;k < STFT_SIZE;k++)
state->OutputAccum[k] += HannWindow[k] * state->FFTbuffer[k].Real /
for(ALsizei k{0};k < STFT_SIZE;k++)
state->OutputAccum[k] += HannWindow[k] * state->FFTbuffer[k].real() /
(0.5 * STFT_HALF_SIZE * OVERSAMP);
/* Shift accumulator, input & output FIFO */
ALsizei j, k;
for(k = 0;k < STFT_STEP;k++) state->OutFIFO[k] = (ALfloat)state->OutputAccum[k];
for(j = 0;k < STFT_SIZE;k++,j++) state->OutputAccum[j] = state->OutputAccum[k];
for(;j < STFT_SIZE;j++) state->OutputAccum[j] = 0.0;
CMakeLists.txt
+1 -1
View File
@@ -755,7 +755,7 @@ ENDIF()
SET(COMMON_OBJS
common/alcomplex.c
common/alcomplex.cpp
common/alcomplex.h
common/align.h
common/almalloc.c
common/alcomplex.c
-122
View File
@@ -1,122 +0,0 @@
#include "config.h"
#include "alcomplex.h"
#include "math_defs.h"
/** Addition of two complex numbers. */
static inline ALcomplex complex_add(ALcomplex a, ALcomplex b)
{
ALcomplex result;
result.Real = a.Real + b.Real;
result.Imag = a.Imag + b.Imag;
return result;
}
/** Subtraction of two complex numbers. */
static inline ALcomplex complex_sub(ALcomplex a, ALcomplex b)
{
ALcomplex result;
result.Real = a.Real - b.Real;
result.Imag = a.Imag - b.Imag;
return result;
}
/** Multiplication of two complex numbers. */
static inline ALcomplex complex_mult(ALcomplex a, ALcomplex b)
{
ALcomplex result;
result.Real = a.Real*b.Real - a.Imag*b.Imag;
result.Imag = a.Imag*b.Real + a.Real*b.Imag;
return result;
}
void complex_fft(ALcomplex *FFTBuffer, ALsizei FFTSize, ALdouble Sign)
{
ALsizei i, j, k, mask, step, step2;
ALcomplex temp, u, w;
ALdouble arg;
/* Bit-reversal permutation applied to a sequence of FFTSize items */
for(i = 1;i < FFTSize-1;i++)
{
for(mask = 0x1, j = 0;mask < FFTSize;mask <<= 1)
{
if((i&mask) != 0)
j++;
j <<= 1;
}
j >>= 1;
if(i < j)
{
temp = FFTBuffer[i];
FFTBuffer[i] = FFTBuffer[j];
FFTBuffer[j] = temp;
}
}
/* Iterative form of Danielson–Lanczos lemma */
for(i = 1, step = 2;i < FFTSize;i<<=1, step<<=1)
{
step2 = step >> 1;
arg = M_PI / step2;
w.Real = cos(arg);
w.Imag = sin(arg) * Sign;
u.Real = 1.0;
u.Imag = 0.0;
for(j = 0;j < step2;j++)
{
for(k = j;k < FFTSize;k+=step)
{
temp = complex_mult(FFTBuffer[k+step2], u);
FFTBuffer[k+step2] = complex_sub(FFTBuffer[k], temp);
FFTBuffer[k] = complex_add(FFTBuffer[k], temp);
}
u = complex_mult(u, w);
}
}
}
void complex_hilbert(ALcomplex *Buffer, ALsizei size)
{
const ALdouble inverse_size = 1.0/(ALdouble)size;
ALsizei todo, i;
for(i = 0;i < size;i++)
Buffer[i].Imag = 0.0;
complex_fft(Buffer, size, 1.0);
todo = size >> 1;
Buffer[0].Real *= inverse_size;
Buffer[0].Imag *= inverse_size;
for(i = 1;i < todo;i++)
{
Buffer[i].Real *= 2.0*inverse_size;
Buffer[i].Imag *= 2.0*inverse_size;
}
Buffer[i].Real *= inverse_size;
Buffer[i].Imag *= inverse_size;
i++;
for(;i < size;i++)
{
Buffer[i].Real = 0.0;
Buffer[i].Imag = 0.0;
}
complex_fft(Buffer, size, -1.0);
}
common/alcomplex.cpp
+76
View File
@@ -0,0 +1,76 @@
#include "config.h"
#include "alcomplex.h"
#include <cmath>
namespace {
constexpr double Pi{3.141592653589793238462643383279502884};
} // namespace
void complex_fft(std::complex<double> *FFTBuffer, int FFTSize, double Sign)
{
/* Bit-reversal permutation applied to a sequence of FFTSize items */
for(int i{1};i < FFTSize-1;i++)
{
int j{0};
for(int mask{1};mask < FFTSize;mask <<= 1)
{
if((i&mask) != 0)
j++;
j <<= 1;
}
j >>= 1;
if(i < j)
std::swap(FFTBuffer[i], FFTBuffer[j]);
}
/* Iterative form of Danielson–Lanczos lemma */
int step{2};
for(int i{1};i < FFTSize;i<<=1, step<<=1)
{
int step2{step >> 1};
double arg{Pi / step2};
std::complex<double> w{std::cos(arg), std::sin(arg)*Sign};
std::complex<double> u{1.0, 0.0};
for(int j{0};j < step2;j++)
{
for(int k{j};k < FFTSize;k+=step)
{
std::complex<double> temp{FFTBuffer[k+step2] * u};
FFTBuffer[k+step2] = FFTBuffer[k] - temp;
FFTBuffer[k] += temp;
}
u *= w;
}
}
}
void complex_hilbert(std::complex<double> *Buffer, int size)
{
const double inverse_size = 1.0/(double)size;
for(int i{0};i < size;i++)
Buffer[i].imag(0.0);
complex_fft(Buffer, size, 1.0);
int todo{size>>1};
int i{0};
Buffer[i++] *= inverse_size;
while(i < todo)
Buffer[i++] *= 2.0*inverse_size;
Buffer[i++] *= inverse_size;
for(;i < size;i++)
Buffer[i] = std::complex<double>{};
complex_fft(Buffer, size, -1.0);
}
common/alcomplex.h
+3 -17
View File
@@ -1,17 +1,7 @@
#ifndef ALCOMPLEX_H
#define ALCOMPLEX_H
#include "AL/al.h"
#ifdef __cplusplus
extern "C" {
#endif
typedef struct ALcomplex {
ALdouble Real;
ALdouble Imag;
} ALcomplex;
#include <complex>
/**
* Iterative implementation of 2-radix FFT (In-place algorithm). Sign = -1 is
@@ -20,7 +10,7 @@ typedef struct ALcomplex {
* FFTBuffer[0...FFTSize-1]. FFTBuffer is an array of complex numbers, FFTSize
* MUST BE power of two.
*/
void complex_fft(ALcomplex *FFTBuffer, ALsizei FFTSize, ALdouble Sign);
void complex_fft(std::complex<double> *FFTBuffer, int FFTSize, double Sign);
/**
* Calculate the complex helical sequence (discrete-time analytical signal) of
@@ -29,10 +19,6 @@ void complex_fft(ALcomplex *FFTBuffer, ALsizei FFTSize, ALdouble Sign);
* Buffer[0...size-1]. Buffer is an array of complex numbers, size MUST BE
* power of two.
*/
void complex_hilbert(ALcomplex *Buffer, ALsizei size);
#ifdef __cplusplus
} // extern "C"
#endif
void complex_hilbert(std::complex<double> *Buffer, int size);
#endif /* ALCOMPLEX_H */