I implemented FFT in C

Question

I wrote Cooley-Tukey's FFT algorithm in C. And I want to know how safe is this code? Is it just trash? How can I make it more memory safe?

#include <stdio.h>
#include <complex.h>
#include <stdlib.h>
#include <math.h>
#define PI 3.14159265358979323846
void fft(double complex *x, unsigned int N, double complex *result){
 if(N == 1){
 result[0] = x[0];
 return;
 }
 double complex sub_even[N/2];
 double complex sub_odd[N/2];
 for(unsigned int i = 0; i < N/2; i++){
 sub_even[i] = x[2 * i];
 sub_odd[i] = x[2 * i + 1];
 }
 
 double complex *even = (double complex *)malloc(sizeof(double complex) * (N/2));
 double complex *odd = (double complex *)malloc(sizeof(double complex) * (N/2));
 fft(sub_even, N/2, even);
 fft(sub_odd, N/2, odd);
 for(unsigned int k = 0; k < N / 2; k++){
 double complex twiddle = cexp(-2 * I * PI * k / N) * odd[k];
 result[k] = even[k] + twiddle;
 result[k + N/2] = even[k] - twiddle;
 }
 free(even);
 free(odd);
 
 
}
int main(){
 unsigned int N = 8;
 double complex x[N];
 for(unsigned int i = 0; i < N; i++){
 x[i] = i+1;
 }
 double complex *result = (double complex *)malloc(sizeof(double complex) * N);
 fft(x, N, result);
 for(unsigned int i = 0; i < N; i++){
 printf("%.2f + %.8fj\n", creal(result[i]), cimag(result[i]));
 }
 return 0;
}

• Jim J. Jewett: What should happen if N is not a power of 2?

Answer 1

double complex sub_even[N/2];
double complex sub_odd[N/2];

These are VLAs and they have the usual problem: it's easy to cause a stack overflow this way. No big deal if N=8, but that's not really a typical FFT size.

These arrays are also not necessary, they only change how the data is indexed (with increasingly large gaps) which you can do with arithmetic instead of by physically moving the data.

double complex *even = (double complex *)malloc(sizeof(double complex) * (N/2));
double complex *odd = (double complex *)malloc(sizeof(double complex) * (N/2));

Even more memory, but FFT can be implemented in-place (transforming the input into the output), or using one buffer of size N if you prefer it that way (just one, allocated before the recursion or pseudo-recursion, not a size N allocation at every level of the recursion).

As for safety, what happens if malloc fails? There's no provision for that (but it would even better to not allocate in the first place). What happens if some huge N is passed in such that sizeof(double complex) * (N/2) = 0 or in general wraps to a small number? Can't happen here really, the VLAs would explode first. But once you remove them, it becomes possible - although I think low-relevance because the only way it happens is if N is wrong (the input array cannot be sufficiently large, its size would also have wrapped), whether intentionally or due to a mistake upstream, either way this code is neither responsible for the problem nor able to do anything about it. But with all that said, it may be a better practice to at least not add an additional bad consequence.

cexp(-2 * I * PI * k / N)

Twiddles can be computed significantly cheaper than exp-ing them one by one.

These basic techniques along with transforming the recursion into iteration are shown for example on wikipedia, although it does not show an efficient way to implement the bit-reverse permutation, for which I refer you to Practically efficient methods for performing bit-reversed permutation in C++11 on the x86-64 architecture (you can do it in C too, and on different architectures).

• RudraSama: Yeah, about those sub_even and sub_odd array, I understood it after running N for 2^18. It gave stack segmentation error LOL. And about in-place (bit reverse) FFT, Yeah I am aware of it, I just wanted to implement it in recursion style. But Thanks for your suggestions, it really helped me to learn alott. And yeah I will try to transform it into iteration. Thanks again!

Answer 2

We assume that N at each level is exactly divisible by 2, but we never confirm this. That can lead to users wondering why up to half of their input is ignored.

---

Consider this allocation:

double complex *result = (double complex *)malloc(sizeof(double complex) * N);

This is C, so the void* returned by malloc() can be assigned to result without any cast. Remove the pointless clutter.

Also, it's a good habit to to use sizeof with an expression rather than a type (*result instead of (double complex)) since that can avoid tediously looking up the target's type when it's far removed from the assignment:

double complex *result = malloc(sizeof *result * N);

And don't forget to ensure that result isn't a null pointer before using it (other than as argument to free()).

---

We're using much more memory than we need to, by allocating at every step of the recursion. We have this lovely result array that we don't use until after the recursive call returns, so that can be used as scratch space. Or re-work to allocate a scratch space at top level and pass that through the recursive calls.

Answer 3

Memory

Seem strange to use VLAs and multiple malloc() allocations - and not checking for success - that is unsafe.

Consider instead 1 allocation and add a check:

// Allocate to the size of the reference object, not the type.
// It is easier to review and maintain.
double complex *dc = malloc(sizeof dc[0] * N/2 * 4);
if (dc == NULL) {
 handle_OOM(); // Out of memory
 ...
} else {
 double complex *sub_even = &dc[0];
 double complex *sub_odd = &dc[N/2];
 ...
 double complex *even = &dc[N/2 * 2];
 double complex *odd = &dc[N/2 * 3];
 ...
 free(dc);
}

Code should check if N is not a power of 2.

Pedantic code would check if N is too large.

if (N > SIZE_MAX/2/sizeof dc[0]) {
 Handle_bad_size();
}

Use const and restrict

When able, use const. It makes for safer maintenance, wider applicability and potential optimizations.

When able, use restrict. Compiler can assume refenced data does not overlap and so allows for potential optimizations. It also lets the user know data should not overlap. As is, x and result could overlap and code works correctly. Yet I suspect later code re-org may not support that. Do you really want the user to be allowed to provide overlapped data?

// void fft(double complex *x, unsigned int N, double complex *result){
void fft(const double complex * restrict x, unsigned int N, double complex *restrict result){

Printing

Avoid %f (fixed point) for general floating point printing. It it verbose when values are large and uninformative when values are small. Consider %g or even %a.

IMO, this improves output as small values do not report like a zero value.

//printf("%.2f + %.8fj\n", creal(result[i]), cimag(result[i]));
printf("%g + %gj\n", creal(result[i]), cimag(result[i]));

Formatting

Save time, use an auto formatter for code.

Sizing

fft() uses unsigned. size_t is the Goldilocks type for array sizing and indexing. size_t is neither too narrow nor too wide for a general purpose sizing with fft().

I'd expect fft() calls to pass in a size_t for size. unsigned may narrow a size_t and thus create a pesky compiler warning.

Design

Perhaps allocate the result space?

// void fft(double complex *x, unsigned int N, double complex *result){
double complex *fft_alloc(size_t n, const double complex x[n]) {

• RudraSama: "Do you really want the user to be allowed to provide overlapped data?" Well it is for my own Mini project, hehe. Btw thanks for soooo much detailed explanation, It will really help me to write safe code next time ughh.... And btw you are right, I should start using code formatter lol. I write very ugly code. Ugh...
• Martin R: I agree about the %g format recommendation, but I wonder what the advantage of using %a would be, unless you want the output to be machine-readable. I don't think that it is "more informative" than %f for humans.
• chux: @MartinR I find %a useful in understanding the least significant bits, especially in comparing a computed answer versus an expected one, the ULP - how close was code? Often I use both %.*g and %a.
• BЈовић: what is an auto formatter?
• chux: @BЈовић An auto formatter is a feature in software applications that automatically formats your text or data, saving you time and ensuring consistency in its appearance, common with many IDEs.

Answer 4

Code reformatted without semantic (or syntactic) changes for better readability:

#include <complex.h>
#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#define PI 3.14159265358979323846
void fft(double complex *x, unsigned int N, double complex *result) {
 if (N == 1) {
 result[0] = x[0];
 return;
 }
 double complex sub_even[N / 2];
 double complex sub_odd[N / 2];
 for (unsigned int i = 0; i < N / 2; i++) {
 sub_even[i] = x[2 * i];
 sub_odd[i] = x[2 * i + 1];
 }
 double complex *even = (double complex *)malloc(sizeof(double complex) * (N / 2));
 double complex *odd = (double complex *)malloc(sizeof(double complex) * (N / 2));
 fft(sub_even, N / 2, even);
 fft(sub_odd, N / 2, odd);
 for (unsigned int k = 0; k < N / 2; k++) {
 double complex twiddle = cexp(-2 * I * PI * k / N) * odd[k];
 result[k] = even[k] + twiddle;
 result[k + N / 2] = even[k] - twiddle;
 }
 free(even);
 free(odd);
}
int main() {
 unsigned int N = 8;
 double complex x[N];
 for (unsigned int i = 0; i < N; i++) {
 x[i] = i + 1;
 }
 double complex *result = (double complex *)malloc(sizeof(double complex) * N);
 fft(x, N, result);
 for (unsigned int i = 0; i < N; i++) {
 printf("%.2f + %.8fj\n", creal(result[i]), cimag(result[i]));
 }
 return 0;
}

Answer 5

You are defining PI and also importing math.h.

M_PI is already defined in math.h, so no need to define it again.

• Andrew Henle: The (draft) C24 standard does not specify any M_PI definition.
• Oddthinking: @AndrewHenle: You are right. It is part of the POSIX standard but not the ISO one.