I'm trying to implement A Radix-5,Radix-3 FFT in C , I already managed to write a Radix-2 but I have some type of bug when it comes to Radix 3 or 5, let's say I do an FFT of 3 samples, that would show the correct results, however if I do FFT of 9 which is 3 * 3, it doesn't show the correct results.
I originally took the code from Python, it works there and I'm trying to simply "copy" it to C , Here is the original Python code:
def fft(x):
"""
radix-2,3,5 FFT algorithm
"""
N = len(x)
if N <= 1:
return x
elif N % 2 == 0:
# For multiples of 2 this formula works
even = fft(x[0::2])
odd = fft(x[1::2])
T = [np.exp(-2j*np.pi*k/N)*odd[k] for k in range(N//2)]
return [even[k] T[k] for k in range(N//2)] \
[even[k] - T[k] for k in range(N//2)]
elif N % 3 == 0:
# Optional, implementing factor 3 decimation
p0 = fft(x[0::3])
p1 = fft(x[1::3])
p2 = fft(x[2::3])
# This will construct the output output without the simplifications
# you can do explorint symmetry
for k in range(N):
return [p0[k % (N//3)]
p1[k % (N//3)] * np.exp(-2j*np.pi*k/N)
p2[k % (N//3)] * np.exp(-4j*np.pi*k/N)]
elif N % 5 == 0:
#factor 5 decimation
p0 = fft(x[0::5])
p1 = fft(x[1::5])
p2 = fft(x[2::5])
p3 = fft(x[3::5])
p4 = fft(x[4::5])
return [p0[k % (N//5)]
p1[k % (N//5)] * np.exp(-2j*np.pi*k/N)
p2[k % (N//5)] * np.exp(-4j*np.pi*k/N)
p3[k % (N//5)] * np.exp(-6j*np.pi*k/N)
p4[k % (N//5)] * np.exp(-8j*np.pi*k/N)
for k in range(N)]
x = [1,1.00071,1.00135,1.00193,1.00245,1.0029,1.00329,1.00361,1.00387]
assert(np.allclose(fft(x), np.fft.fft(x)))
And here is my code in C (fft.hpp):
#define _USE_MATH_DEFINES
#pragma once
#include <cmath>
#include <vector>
#include <complex>
using std::vector;
using std::complex;
vector<complex<float>> slicing(vector<complex<float>> vec, unsigned int X, unsigned int Y, unsigned int stride)
{
// To store the sliced vector
vector<complex<float>> result;
// Copy vector using copy function()
int i = X;
while (result.size() < Y)
{
result.push_back(vec[i]);
i = i stride;
}
// Return the final sliced vector
return result;
}
void fft(vector<complex<float>>& x)
{
// Check if it is splitted enough
const size_t N = x.size();
if (N <= 1)
return;
else if (N % 2 == 0)
{
//Radix-2
vector<complex<float>> even = slicing(x, 0, N / 2, 2); //split the inputs in even / odd indices subarrays
vector<complex<float>> odd = slicing(x, 1, N / 2, 2);
// conquer
fft(even);
fft(odd);
// combine
for (size_t k = 0; k < N / 2; k)
{
complex<float> t = std::polar<float>(1.0, -2 * M_PI * k / N) * odd[k];
x[k] = even[k] t;
x[k N / 2] = even[k] - t;
}
}
else if (N % 3 == 0)
{
//Radix-3
//factor 3 decimation
vector<complex<float>> p0 = slicing(x, 0, N / 3, 3);
vector<complex<float>> p1 = slicing(x, 1, N / 3, 3);
vector<complex<float>> p2 = slicing(x, 2, N / 3, 3);
fft(p0);
fft(p1);
fft(p2);
for (int i = 0; i < N; i )
{
complex<float> temp = p0[i % (int)N / 3];
temp = (p1[i % (int)N / 3] * std::polar<float>(1.0, -2 * M_PI * i / N));
temp = (p2[i % (int)N / 3] * std::polar<float>(1.0, -4 * M_PI * i / N));
x[i] = temp;
}
}
else if (N % 5 == 0)
{
//Radix-5
//factor 5 decimation
vector<complex<float>> p0 = slicing(x, 0, N / 5, 5);
vector<complex<float>> p1 = slicing(x, 1, N / 5, 5);
vector<complex<float>> p2 = slicing(x, 2, N / 5, 5);
vector<complex<float>> p3 = slicing(x, 3, N / 5, 5);
vector<complex<float>> p4 = slicing(x, 4, N / 5, 5);
fft(p0);
fft(p1);
fft(p2);
fft(p3);
fft(p4);
for (int i = 0; i < N; i )
{
complex<float> temp = p0[i % (int)N / 5];
temp = (p1[i % (int)N / 5] * std::polar<float>(1.0, -2 * M_PI * i / N));
temp = (p2[i % (int)N / 5] * std::polar<float>(1.0, -4 * M_PI * i / N));
temp = (p3[i % (int)N / 5] * std::polar<float>(1.0, -6 * M_PI * i / N));
temp = (p4[i % (int)N / 5] * std::polar<float>(1.0, -8 * M_PI * i / N));
x[i] = temp;
}
}
}
And main.cpp:
#define _USE_MATH_DEFINES
#include <stdio.h>
#include <iostream>
#include "fft.hpp"
typedef vector<complex<float>> complexSignal;
int main()
{
complexSignal abit;
int N = 9;
abit.push_back({1,0});
abit.push_back({1.00071 ,0 });
abit.push_back({1.00135 ,0 });
abit.push_back({1.00193 ,0 });
abit.push_back({1.00245 ,0 });
abit.push_back({1.0029 ,0 });
abit.push_back({1.00329 ,0 });
abit.push_back({1.00361 ,0 });
abit.push_back({1.00387 ,0 });
//abit.push_back({ 1.0029 ,0 });
std::cout << "Before:" << std::endl;
for (int i = 0; i < N; i )
{
//abit.push_back(data[0][i]);
std::cout << abit[i] << std::endl;
}
std::cout << "After:" << std::endl;
fft(abit);
for (int i = 0; i < N; i )
{
std::cout << abit[i] << std::endl;
}
return 0;
}
I'm getting the following output from CPP:
(9.02011,0)
(5.83089,-4.89513)
(0.700632,-3.98993)
(-0.000289979,0.000502368)
(-0.00218513,0.000362784)
(-0.00179241,0.00139188)
(-0.000289979,-0.000502368)
(0.000175771,-0.00354373)
(-0.003268,-0.00558837)
While I should get:
(9.020109999999999 0j)
(-0.0032675770104925446 0.005588577982060319j)
(-0.0023772289746976797 0.0024179090499282354j)
(-0.0022250000000012538 0.0011691342951078987j)
(-0.002185194014811494 0.00036271471530890747j)
(-0.0021851940148113033-0.00036271471530980844j)
(-0.0022249999999994774-0.0011691342951105632j)
(-0.002377228974696629-0.0024179090499291786j)
(-0.00326757701049002-0.005588577982061138j)
As you can see only the first result is correct.
CodePudding user response：
In Python it says k % (N//3), in C it says i % (int)N / 3. These are not the same! See operator precedence in C . You need to use parenthesis to execute the operations in the right order: i % ((int)N / 3).
Pepijn Kramer’s comment about changing
-2 * M_PI * i / N
into
-2.0 * M_PI * static_cast<double>(i)/static_cast<double>(N)
is good practice, in my opinion, but shouldn’t make a difference in this case because of the left-to-right evaluation and the automatic promotion of integers to floats in a float context. Still, ensuring all your constants are floats where needed doesn’t cost much effort and can prevent difficult to spot bugs.