wip : fftw without std::complex

This commit is contained in:
Georgi Gerganov
2020-12-05 14:00:32 +02:00
parent 72d0ca630d
commit cf35ed33c9
2 changed files with 48 additions and 30 deletions

View File

@@ -28,31 +28,48 @@ int reverse(int N, int n) {
return p;
}
void ordina(std::complex<float>* f1, int N) {
std::complex<float> f2[GGWave::kMaxSamplesPerFrame];
for(int i = 0; i < N; i++)
f2[i] = f1[reverse(N, i)];
for(int j = 0; j < N; j++)
f1[j] = f2[j];
void ordina(float * f1, int N) {
float f2[2*GGWave::kMaxSamplesPerFrame];
for (int i = 0; i < N; i++) {
int ir = reverse(N, i);
f2[2*i + 0] = f1[2*ir + 0];
f2[2*i + 1] = f1[2*ir + 1];
}
for (int j = 0; j < N; j++) {
f1[2*j + 0] = f2[2*j + 0];
f1[2*j + 1] = f2[2*j + 1];
}
}
void transform(std::complex<float>* f, int N) {
void transform(float * f, int N) {
ordina(f, N); //first: reverse order
std::complex<float> *W;
W = (std::complex<float> *)malloc(N / 2 * sizeof(std::complex<float>));
W[1] = std::polar(1., -2. * M_PI / N);
W[0] = 1;
for(int i = 2; i < N / 2; i++)
W[i] = pow(W[1], i);
float * W;
W = (float *)malloc(N*sizeof(float));
W[2*1 + 0] = cos(-2.*M_PI/N);
W[2*1 + 1] = sin(-2.*M_PI/N);
W[2*0 + 0] = 1;
W[2*0 + 1] = 0;
for (int i = 2; i < N / 2; i++) {
W[2*i + 0] = cos(-2.*i*M_PI/N);
W[2*i + 1] = sin(-2.*i*M_PI/N);
}
int n = 1;
int a = N / 2;
for(int j = 0; j < log2(N); j++) {
for(int i = 0; i < N; i++) {
if(!(i & n)) {
std::complex<float> temp = f[i];
std::complex<float> Temp = W[(i * a) % (n * a)] * f[i + n];
f[i] = temp + Temp;
f[i + n] = temp - Temp;
int wi = (i * a) % (n * a);
int fi = i + n;
float a = W[2*wi + 0];
float b = W[2*wi + 1];
float c = f[2*fi + 0];
float d = f[2*fi + 1];
float temp[2] = { f[2*i + 0], f[2*i + 1] };
float Temp[2] = { a*c - b*d, b*c + a*d };
f[2*i + 0] = temp[0] + Temp[0];
f[2*i + 1] = temp[1] + Temp[1];
f[2*fi + 0] = temp[0] - Temp[0];
f[2*fi + 1] = temp[1] - Temp[1];
}
}
n *= 2;
@@ -61,16 +78,18 @@ void transform(std::complex<float>* f, int N) {
free(W);
}
void FFT(std::complex<float>* f, int N, float d) {
void FFT(float * f, int N, float d) {
transform(f, N);
for(int i = 0; i < N; i++)
f[i] *= d; //multiplying by step
for (int i = 0; i < N; i++) {
f[2*i + 0] *= d;
f[2*i + 1] *= d;
}
}
void FFT(float * src, std::complex<float>* dst, int N, float d) {
void FFT(float * src, float * dst, int N, float d) {
for (int i = 0; i < N; ++i) {
dst[i].real(src[i]);
dst[i].imag(0);
dst[2*i + 0] = src[i];
dst[2*i + 1] = 0.0f;
}
FFT(dst, N, d);
}
@@ -177,8 +196,8 @@ bool GGWave::init(int textLength, const char * stext, const TxProtocol & aProtoc
m_rxData.fill(0);
for (int i = 0; i < m_samplesPerFrame; ++i) {
m_fftOut[i].real(0.0f);
m_fftOut[i].imag(0.0f);
m_fftOut[2*i + 0] = 0.0f;
m_fftOut[2*i + 1] = 0.0f;
}
return true;
@@ -378,7 +397,7 @@ void GGWave::receive(const CBDequeueAudio & CBDequeueAudio) {
double fsum = 0.0;
for (int i = 0; i < m_samplesPerFrame; ++i) {
m_sampleSpectrum[i] = (m_fftOut[i].real()*m_fftOut[i].real() + m_fftOut[i].imag()*m_fftOut[i].imag());
m_sampleSpectrum[i] = (m_fftOut[2*i + 0]*m_fftOut[2*i + 0] + m_fftOut[2*i + 1]*m_fftOut[2*i + 1]);
fsum += m_sampleSpectrum[i];
}
for (int i = 1; i < m_samplesPerFrame/2; ++i) {
@@ -440,7 +459,7 @@ void GGWave::receive(const CBDequeueAudio & CBDequeueAudio) {
FFT(m_fftIn.data(), m_fftOut.data(), m_samplesPerFrame, 1.0);
for (int i = 0; i < m_samplesPerFrame; ++i) {
m_sampleSpectrum[i] = (m_fftOut[i].real()*m_fftOut[i].real() + m_fftOut[i].imag()*m_fftOut[i].imag());
m_sampleSpectrum[i] = (m_fftOut[2*i + 0]*m_fftOut[2*i + 0] + m_fftOut[2*i + 1]*m_fftOut[2*i + 1]);
}
for (int i = 1; i < m_samplesPerFrame/2; ++i) {
m_sampleSpectrum[i] += m_sampleSpectrum[m_samplesPerFrame - i];