buddhabrot
/
fusion-zauberstab
1
0
Fork 1
Code Issues Pull Requests Projects Releases Wiki Activity

refactor fft to be type generic 

Signed-off-by: Thomas Schmid <tom@lfence.de>
This commit is contained in:
Thomas 2022-06-19 14:24:02 +02:00
parent 5140512b56
commit 3d8e9deddd
3 changed files with 80 additions and 80 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

74
firmware/include/fft.h
View File

@ -2,5 +2,75 @@
#include <complex>
void fft(std::complex<float> *input, std::complex<float> *output, uint32_t N);
void rfft(std::complex<float> *input, std::complex<float> *output, uint32_t N);
template <class T>
struct FFT
{
static void fft(std::complex<T> *samples, std::complex<T> *output, uint32_t N)
{
uint8_t log2n = (uint8_t)std::log2(N) + 0.5f;
std::complex<T> I(0.0, 1.0);
if (N == 1)
{
output[0] = samples[0];
return;
}
// shuffle array
for (int i = 0; i < N; i++)
{
output[i] = samples[FFT::bitReverse(i, log2n)];
}
for (int s = 1; s <= log2n; s++)
{
uint32_t m = 1 << s; // 2^s
std::complex<T> wm = std::exp(-2.0f * (T)M_PI * I / (std::complex<T>)m);
for (int k = 0; k < N; k += m)
{
std::complex<T> w = 1.f;
for (int j = 0; j < m / 2; j++)
{
std::complex<T> t = w * output[k + j + m / 2];
std::complex<T> u = output[k + j];
output[k + j] = u + t;
output[k + j + m / 2] = u - t;
w = w * wm;
}
}
}
}
static void rfft(std::complex<T> *input, std::complex<T> *output, uint32_t N)
{
std::complex<T> I(0.0, 1.0);
for (int i = 0; i < N / 2; i++)
{
input[i] = input[i] + I * input[i + N / 2];
}
FFT<T>::fft(input, output, N / 2);
for (int i = 0; i < N / 2; i++)
{
output[i] = (output[i] + std::conj(output[(N / 2) - i])) / 2.;
}
for (int i = N / 2; i < N; i++)
{
output[i] = -I * (output[i] - std::conj(output[(N / 2) - i])) / 2.;
}
}
private:
static unsigned int bitReverse(unsigned int x, int log2n)
{
int n = 0;
for (int i = 0; i < log2n; i++)
{
n <<= 1;
n |= (x & 1);
x >>= 1;
}
return n;
}
};

16
firmware/src/applications/fft_test.cpp
View File

@ -4,13 +4,13 @@
#define N_SAMPLES 256
std::complex<float> samples[N_SAMPLES];
std::complex<float> z[N_SAMPLES];
double vReal[N_SAMPLES];
double vImag[N_SAMPLES];
uint32_t sample_counter = 0;
unsigned long max_dt = 0;
unsigned long last_sample = 0;
static std::complex<float> samples[N_SAMPLES];
static std::complex<float> z[N_SAMPLES];
static double vReal[N_SAMPLES];
static double vImag[N_SAMPLES];
static uint32_t sample_counter = 0;
static unsigned long max_dt = 0;
static unsigned long last_sample = 0;
void FFTTestApp::init()
{
@ -37,7 +37,7 @@ void FFTTestApp::loop()
if (sample_counter == N_SAMPLES)
{
unsigned long start = micros();
fft(samples, z, N_SAMPLES);
FFT<float>::fft(samples, z, N_SAMPLES);
unsigned long end = micros();
float max = 0.f;

70
firmware/src/fft.cpp
View File

@ -1,70 +0,0 @@
#include "fft.h"
#include "math.h"
unsigned int bitReverse(unsigned int x, int log2n)
{
int n = 0;
for (int i = 0; i < log2n; i++)
{
n <<= 1;
n |= (x & 1);
x >>= 1;
}
return n;
}
void fft(std::complex<float> *samples, std::complex<float> *output, uint32_t N)
{
uint8_t log2n = (uint8_t)std::log2(N) + 0.5f;
std::complex<float> I(0.0, 1.0);
if (N == 1)
{
output[0] = samples[0];
return;
}
// shuffle array
for (int i = 0; i < N; i++)
{
output[i] = samples[bitReverse(i, log2n)];
}
for (int s = 1; s <= log2n; s++)
{
uint32_t m = 1 << s; // 2^s
std::complex<float> wm = std::exp(-2.0f * (float)M_PI * I / (std::complex<float>)m);
for (int k = 0; k < N; k += m)
{
std::complex<float> w = 1.f;
for (int j = 0; j < m / 2; j++)
{
std::complex<float> t = w * output[k + j + m / 2];
std::complex<float> u = output[k + j];
output[k + j] = u + t;
output[k + j + m / 2] = u - t;
w = w * wm;
}
}
}
}
void rfft(std::complex<float> *input, std::complex<float> *output, uint32_t N)
{
std::complex<float> I(0.0, 1.0);
for (int i = 0; i < N / 2; i++)
{
input[i] = input[i] + I * input[i + N / 2];
}
fft(input, output, N / 2);
for (int i = 0; i < N / 2; i++)
{
output[i] = (output[i] + std::conj(output[(N / 2) - i])) / 2.f;
}
for (int i = N / 2; i < N; i++)
{
output[i] = -I * (output[i] - std::conj(output[(N / 2) - i])) / 2.f;
}
}