Cant figure out how to set up re-cursive functions

1

#![allow(non_snake_case)]

use num_complex::Complex64;

fn main() {
    let one= Complex64::new(1.0, 0.0);
    let zz = Complex64::new(0.0, 0.0);
    let mut data = vec![one,one,one,one, zz,zz,zz,zz,];
    let     test = vec![one,one,one,one, zz,zz,zz,zz,];

    fft (&mut data); // ifft (&mut data, DC);

    let mut E = 0_u32;  for (k,s) in data.iter().enumerate() { if (s -test[k]).norm() > 1.0e-8 { E += 1; }  }
    println!(" abs(ifft vs test data) Errors {}", E);
}

const PI: f64 = 3.14159265358979323846264338327950288_f64;

fn fft (x: &mut Vec<Complex64> )
{   let N = x.len() as usize; if N <= 1 { return ; }
    let mut even = vec![Complex64::default(); N/2];
    let mut odd  = vec![Complex64::default(); N/2];
    let mut k = 0_usize;
    while k< N/2 { even.push(x[2*k]); odd.push(x[2*k+1]); k += 1; }
///////////
    assert!( even[0]==x[0] && (odd[0]==x[1])) ;
///////////
    fft( &mut even); fft( &mut odd);
    // for k in 0..N/2
    k= 0; let mut Ninv = 0.0_f64;
    while k< N/2
    {   let t= odd[k]* (Complex64::from_polar(1.0, -2.0*PI*Ninv) ).exp(); // (k as f64)/(N as f64)) ).exp();
        x[k]     = even[k] +t;
        x[k +N/2]= even[k] -t;
        Ninv += 1.0/(N as f64);
        k += 1;
    }
}
2

The code you have posted is incomplete (e.g. doesn't include import/definition of Complex64), and you have not included the error message, so I am not able to see the details of the problem.

Do you get an error message? What is it?

3

the assert fails

4

Looks like you are comparing floating point numbers for equality. One should never do this because the same values calculated by different means using floating point are very likely to not be equal. I would expect your asserts to fail.

This perhaps counter intuitive result is due to the inherent inaccuracies of floating point, accumulation of rounding errors and so on.

See "What Every Computer Scientist Should Know About Floating-Point Arithmetic": https://www.itu.dk/~sestoft/bachelor/IEEE754_article.pdf

5

But he hasn't done any arithmetic on the numbers. He's just copying them to a different array. Perhaps he has a Nan.

6

Ah yes. I should have looked closer.

7

This allocates two vectors and fills them with the default value. Their initial length is N/2.

    let mut even = vec![Complex64::default(); N/2];
    let mut odd  = vec![Complex64::default(); N/2];

Then you push some other values onto the end of the vectors, into positions N/2, N/2+1, etc.

    let mut k = 0_usize;
    while k< N/2 { even.push(x[2*k]); odd.push(x[2*k+1]); k += 1; }

But your assert looks at position 0.

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

This assert succeeds (but the stack overflow for some other reason):

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

Perhaps what you wanted was

let mut even = Vec::with_capacity(N/2);
let mut odd = Vec::with_capacity(N/2);

...but I haven't analyzed the algorithm at all, so that's just a guess.

8

Rather that copying odds and evens into separate arrays wouldn't it be easier, faster, to iterate over the elements with a stride: itertools::Stride itertools - Rust ?

9

They're porting C or C++ code as I recall, but I agree a "Rustification" type review would be in order at some point in the process.

10

After if N <= 1 {return;} we have N > 1 and thus N/2 > 0. Then even[0] = Complex64::default(), but x[0] is arbitrary, in this case one.

11

the floats should be dead equal...even[0] is just a copy of x[0],...ditto for odd[0] and x[1]

12

I presume that code is supposed to be an implementation of the. recursive Cooley-Tukey as shown in the pseudo code here: Cooley–Tukey FFT algorithm - Wikipedia Where the function takes an array, a length and a stride.

Looks like it should be a nice exercise to implement that in Rust using the itertools::stride.

13

That's not true; reread my longer post. even[N/2] is a copy of x[0] and odd[N/2] is a copy of x[1]. I'm pretty sure you want Vec::with_capacity(N/2) and not vec![Complex64::default(); N/2] after running it on the playground.

14

Is there an ownership issue in the while { } ??

15

Vec::with_capacity(N/2)
gets me passed the assert ...for reasons I dont understand...I don`t see how the assert passes ::with_cap.. and not the original stmt
...Also the final while { } has bugs...and I not sure about the re-cursive calls

16

@ZiCog ...You are correct ...it is the C-T FFT algorithm...I converting some of my scientific software to Rust in hopes of building up my Rust savvy

17 
    let mut data = vec![one,one,one,one, zz,zz,zz,zz,];
    let     test = vec![one,one,one,one, zz,zz,zz,zz,];

The test values are not the actual result of FFT on the data values.

let t= odd[k]* (Complex64::from_polar(1.0, -2.0*PI*Ninv) ).exp();

You get the .exp() as part of the conversion from polar (you should remove the .exp()).

Here's what I got after those changes and some non-algorithmic changes. Don't look yet if you'd rather work out the bugs yourself. I didn't test anything except the singular example. I didn't attempt an in-place version.

18

You're comparing the default float to x[0]. The vec! macro allocates and populates N/2 items, and then your code pushes an additional N/2 items.

19

And here's a snippet to illustrate the differences between vec![...], Vec::new(), and Vec::with_capacity().

20

@nbaraz: aaahhhh yes, yes...I get it now