Cant figure out how to set up re-cursive functions

@qinedot: correct! the .exp() should be removed

Don't worry. I looked. I'm never going to live long enough to understand what is written there.

Is it possible to write a version of the recursive Cooley-Tukey that looks like the pseudo code linked above so that it's recognisable/understandable?

X0,...,Nāˆ’1 ā† ditfft2(x, N, s):                         DFT of (x0, xs, x2s, ..., x(N-1)s):
    if N = 1 then
        X0 ā† x0                                        trivial size-1 DFT base case
    else
        X0,...,N/2āˆ’1 ā† ditfft2(x, N/2, 2s)             DFT of (x0, x2s, x4s, ...)
        XN/2,...,Nāˆ’1 ā† ditfft2(x+s, N/2, 2s)           DFT of (xs, xs+2s, xs+4s, ...)
        for k = 0 to N/2āˆ’1 do                          combine DFTs of two halves into full DFT:
            t ā† Xk
            Xk ā† t + exp(āˆ’2Ļ€i k/N) Xk+N/2
            Xk+N/2 ā† t āˆ’ exp(āˆ’2Ļ€i k/N) Xk+N/2
        end for
    end if

Everything's relative I guess. I had to squint at the pseudo-code to follow it, especially the vague list assignment expression, but it is what I used to touch up the original code. Most of my changes were rewriting loops as iterator chains; feel free to ask if you want me to walk through anything and associate it to the pseudo-code.

use num_complex::Complex64 as C64;
use std::f64::consts::PI;
const I: C64 = C64::new(0.0, 1.0);

fn step2(x: &[C64], start: usize) -> Vec<C64> {
    x[start..].iter().step_by(2).cloned().collect()
}

#[allow(non_snake_case)]
fn fft(x: &[C64]) -> Vec<C64> {
   if x.len() <= 1 {
      x.to_vec()
   }  else {
      let (even, odd) = (fft(&step2(x, 0)), fft(&step2(x, 1)));
      let N = x.len();
      let t: Vec<C64> = (0..N/2)
          .map(|k| C64::exp(-2.0*I*PI*(k as f64)/(N as f64))*odd[k])
          .collect();
      (0..N/2).map(|k| even[k] + t[k]).chain(
      (0..N/2).map(|k| even[k] - t[k])).collect()
   }
}


fn main() {
    let data: Vec<C64> = [1, 1, 1, 1, 0, 0, 0, 0].iter()
        .map(|&x| C64::new(f64::from(x), 0.0)).collect();
    println!("{:#?}", fft(&data));
}

Offtopic but something that caught my eye:

assert!( even[0]==x[0] && (odd[0]==x[1])) ;

When using asserts I always recommend using the shortest, most specific assert possible. If you use expressions you will never learn which part of the expression is not true.

If you rewrite the above as

assert_eq!( even[0], x[0] );
assert_eq!( odd[0], x[1] ) ;

it'll be much clearer what and why fails and you'll be in a better position to understand the problem. Copy&pasting the output could also help others help you.

Otherwise assert messages can get very messy, like this one I encountered a couple of years ago:

malloc.c:2369: sysmalloc: Assertion `(old_top == (((mbinptr) (((char *) &((av)->bins[((1) - 1) * 2])) - __builtin_offsetof (struct malloc_chunk, fd)))) && old_size == 0) || ((unsigned long) (old_size) >= (unsigned long)((((__builtin_offsetof (struct malloc_chunk, fd_nextsize))+((2 * (sizeof(size_t))) - 1)) & ~((2 * (sizeof(size_t))) - 1))) && ((old_top)->size & 0x1) && ((unsigned long)old_end & pagemask) == 0)' failed.
4 Likes

@ZiCog: I believe my FFT and inverse FFT code is working now and verified by several test cases (my version of test code). This code is quite compact and provides access to two C-T algorithms.
I'm a Rust newbie, and not ready to get into creation of Rust crates, but I'd be glad to share the code with anyone interested in improving 'newbie code` .