Use a span for the complex_fft/hilbert functions
This commit is contained in:
@@ -171,7 +171,7 @@ void FshifterState::process(const ALsizei samplesToDo, const FloatBufferLine *RE
|
||||
}
|
||||
|
||||
/* Processing signal by Discrete Hilbert Transform (analytical signal). */
|
||||
complex_hilbert(mAnalytic, HIL_SIZE);
|
||||
complex_hilbert(mAnalytic);
|
||||
|
||||
/* Windowing and add to output accumulator */
|
||||
for(k = 0;k < HIL_SIZE;k++)
|
||||
|
||||
@@ -230,7 +230,7 @@ void PshifterState::process(const ALsizei samplesToDo, const FloatBufferLine *RE
|
||||
|
||||
/* ANALYSIS */
|
||||
/* Apply FFT to FFTbuffer data */
|
||||
complex_fft(mFFTbuffer, STFT_SIZE, -1.0);
|
||||
complex_fft(mFFTbuffer, -1.0);
|
||||
|
||||
/* Analyze the obtained data. Since the real FFT is symmetric, only
|
||||
* STFT_HALF_SIZE+1 samples are needed.
|
||||
@@ -302,7 +302,7 @@ void PshifterState::process(const ALsizei samplesToDo, const FloatBufferLine *RE
|
||||
mFFTbuffer[k] = complex_d{};
|
||||
|
||||
/* Apply iFFT to buffer data */
|
||||
complex_fft(mFFTbuffer, STFT_SIZE, 1.0);
|
||||
complex_fft(mFFTbuffer, 1.0);
|
||||
|
||||
/* Windowing and add to output */
|
||||
for(ALsizei k{0};k < STFT_SIZE;k++)
|
||||
|
||||
+29
-29
@@ -4,6 +4,7 @@
|
||||
#include "alcomplex.h"
|
||||
|
||||
#include <cmath>
|
||||
#include <algorithm>
|
||||
|
||||
namespace {
|
||||
|
||||
@@ -11,13 +12,14 @@ constexpr double Pi{3.141592653589793238462643383279502884};
|
||||
|
||||
} // namespace
|
||||
|
||||
void complex_fft(std::complex<double> *FFTBuffer, int FFTSize, double Sign)
|
||||
void complex_fft(const al::span<std::complex<double>> buffer, const double sign)
|
||||
{
|
||||
const size_t fftsize{buffer.size()};
|
||||
/* Bit-reversal permutation applied to a sequence of FFTSize items */
|
||||
for(int i{1};i < FFTSize-1;i++)
|
||||
for(size_t i{1u};i < fftsize-1;i++)
|
||||
{
|
||||
int j{0};
|
||||
for(int mask{1};mask < FFTSize;mask <<= 1)
|
||||
size_t j{0u};
|
||||
for(size_t mask{1u};mask < fftsize;mask <<= 1)
|
||||
{
|
||||
if((i&mask) != 0)
|
||||
j++;
|
||||
@@ -26,25 +28,25 @@ void complex_fft(std::complex<double> *FFTBuffer, int FFTSize, double Sign)
|
||||
j >>= 1;
|
||||
|
||||
if(i < j)
|
||||
std::swap(FFTBuffer[i], FFTBuffer[j]);
|
||||
std::swap(buffer[i], buffer[j]);
|
||||
}
|
||||
|
||||
/* Iterative form of DanielsonLanczos lemma */
|
||||
int step{2};
|
||||
for(int i{1};i < FFTSize;i<<=1, step<<=1)
|
||||
size_t step{2u};
|
||||
for(size_t i{1u};i < fftsize;i<<=1, step<<=1)
|
||||
{
|
||||
int step2{step >> 1};
|
||||
const size_t step2{step >> 1};
|
||||
double arg{Pi / step2};
|
||||
|
||||
std::complex<double> w{std::cos(arg), std::sin(arg)*Sign};
|
||||
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(size_t j{0};j < step2;j++)
|
||||
{
|
||||
for(int k{j};k < FFTSize;k+=step)
|
||||
for(size_t k{j};k < fftsize;k+=step)
|
||||
{
|
||||
std::complex<double> temp{FFTBuffer[k+step2] * u};
|
||||
FFTBuffer[k+step2] = FFTBuffer[k] - temp;
|
||||
FFTBuffer[k] += temp;
|
||||
std::complex<double> temp{buffer[k+step2] * u};
|
||||
buffer[k+step2] = buffer[k] - temp;
|
||||
buffer[k] += temp;
|
||||
}
|
||||
|
||||
u *= w;
|
||||
@@ -52,25 +54,23 @@ void complex_fft(std::complex<double> *FFTBuffer, int FFTSize, double Sign)
|
||||
}
|
||||
}
|
||||
|
||||
void complex_hilbert(std::complex<double> *Buffer, int size)
|
||||
void complex_hilbert(const al::span<std::complex<double>> buffer)
|
||||
{
|
||||
const double inverse_size = 1.0/static_cast<double>(size);
|
||||
std::for_each(buffer.begin(), buffer.end(), [](std::complex<double> &c) { c.imag(0.0); });
|
||||
|
||||
for(int i{0};i < size;i++)
|
||||
Buffer[i].imag(0.0);
|
||||
complex_fft(buffer, 1.0);
|
||||
|
||||
complex_fft(Buffer, size, 1.0);
|
||||
const double inverse_size = 1.0/static_cast<double>(buffer.size());
|
||||
auto bufiter = buffer.begin();
|
||||
const auto halfiter = bufiter + (buffer.size()>>1);
|
||||
|
||||
int todo{size>>1};
|
||||
int i{0};
|
||||
*bufiter *= inverse_size; ++bufiter;
|
||||
bufiter = std::transform(bufiter, halfiter, bufiter,
|
||||
[inverse_size](const std::complex<double> &c) -> std::complex<double>
|
||||
{ return c * (2.0*inverse_size); });
|
||||
*bufiter *= inverse_size; ++bufiter;
|
||||
|
||||
Buffer[i++] *= inverse_size;
|
||||
while(i < todo)
|
||||
Buffer[i++] *= 2.0*inverse_size;
|
||||
Buffer[i++] *= inverse_size;
|
||||
std::fill(bufiter, buffer.end(), std::complex<double>{});
|
||||
|
||||
for(;i < size;i++)
|
||||
Buffer[i] = std::complex<double>{};
|
||||
|
||||
complex_fft(Buffer, size, -1.0);
|
||||
complex_fft(buffer, -1.0);
|
||||
}
|
||||
|
||||
+9
-9
@@ -3,22 +3,22 @@
|
||||
|
||||
#include <complex>
|
||||
|
||||
#include "alspan.h"
|
||||
|
||||
/**
|
||||
* Iterative implementation of 2-radix FFT (In-place algorithm). Sign = -1 is
|
||||
* FFT and 1 is iFFT (inverse). Fills FFTBuffer[0...FFTSize-1] with the
|
||||
* Discrete Fourier Transform (DFT) of the time domain data stored in
|
||||
* FFTBuffer[0...FFTSize-1]. FFTBuffer is an array of complex numbers, FFTSize
|
||||
* MUST BE power of two.
|
||||
* FFT and 1 is iFFT (inverse). Fills the buffer with the Discrete Fourier
|
||||
* Transform (DFT) of the time domain data stored in the buffer. The buffer is
|
||||
* an array of complex numbers, and MUST BE power of two.
|
||||
*/
|
||||
void complex_fft(std::complex<double> *FFTBuffer, int FFTSize, double Sign);
|
||||
void complex_fft(const al::span<std::complex<double>> buffer, const double sign);
|
||||
|
||||
/**
|
||||
* Calculate the complex helical sequence (discrete-time analytical signal) of
|
||||
* the given input using the discrete Hilbert transform (In-place algorithm).
|
||||
* Fills Buffer[0...size-1] with the discrete-time analytical signal stored in
|
||||
* Buffer[0...size-1]. Buffer is an array of complex numbers, size MUST BE
|
||||
* power of two.
|
||||
* Fills the buffer with the discrete-time analytical signal stored in the
|
||||
* buffer. The buffer is an array of complex numbers and MUST BE power of two.
|
||||
*/
|
||||
void complex_hilbert(std::complex<double> *Buffer, int size);
|
||||
void complex_hilbert(const al::span<std::complex<double>> buffer);
|
||||
|
||||
#endif /* ALCOMPLEX_H */
|
||||
|
||||
Reference in New Issue
Block a user