rft 0.1.0

FFT handling library.
Documentation 
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128

use num;
use strided::{Stride, MutStride};
use std::f64::consts::PI;

use {Precision, Complex, ComplexMut};

fn fft<CI: Complex, CO: ComplexMut>(direction: Precision, input: Stride<CI>, mut output: MutStride<CO>) {
	// cache the length
	let length = input.len();

	// base case: the DFT of a single element is itself.
	if length == 1 {
		output[0].set(&input[0]);

		return;
	}

	// split the input into two arrays of alternating elements ("decimate in
	// time")
	let (evens, odds) = input.substrides2();

	// break the output into two halves (front and back, not alternating)
	let (mut left, mut right) = output.split_at_mut(length >> 1);

	// recursively perform two FFTs on alternating elements of the input, writing
	// the results into the first and second half of the output array
	// respectively
	fft(direction, evens, left.reborrow());
	fft(direction, odds, right.reborrow());

	// exp(-2πi/N)
	let twiddle = num::Complex::from_polar(&1.0,
		&(direction * PI as Precision / length as Precision));

	let mut factor = num::Complex::new(1.0, 0.0);

	// combine the subFFTs with the relations:
	//   X_k       = E_k + exp(-2πki/N) * O_k
	//   X_{k+N/2} = E_k - exp(-2πki/N) * O_k
	for (even, odd) in left.iter_mut().zip(right.iter_mut()) {
		let twiddled = factor * odd.to_num();
		let e        = even.to_num();

		even.set(&(e + twiddled));
		odd.set(&(e - twiddled));

		factor = factor * twiddle;
	}
}

/// Applies a forward Cooley-Tukey Fourier transform on the given input and
/// puts the result in the given output.
#[inline(always)]
pub fn forward<CI: Complex, CO: ComplexMut>(input: Stride<CI>, output: MutStride<CO>) {
	// input and output buffers need to be the same length
	debug_assert_eq!(input.len(), output.len());

	// the length has to be a power of two
	debug_assert!(input.len().is_power_of_two(), "length is not a power of two");

	fft(-2.0, input, output);
}

/// Applies an inverse Cooley-Tukey Fourier transform on the given input and
/// puts the result in the given output.
///
/// Note the result is not scaled.
#[inline(always)]
pub fn inverse<CI: Complex, CO: ComplexMut>(input: Stride<CI>, output: MutStride<CO>) {
	// input and output buffers need to be the same length
	debug_assert_eq!(input.len(), output.len());

	// the length has to be a power of two
	debug_assert!(input.len().is_power_of_two(), "length is not a power of two");

	fft(2.0, input, output);
}

#[cfg(test)]
mod tests {
	use num::Complex;
	use strided::{Stride, MutStrided};
	use ::ComplexMut;

	macro_rules! fix {
		($a:expr) => (
			if $a == "-0.00" {
				"0.00".to_owned()
			}
			else {
				$a
			}
		)
	}

	macro_rules! assert_approx_eq {
		($a:expr, $b:expr) => (
			assert_eq!(fix!(format!("{:.2}", $a.re)), fix!(format!("{:.2}", $b.re)));
			assert_eq!(fix!(format!("{:.2}", $a.im)), fix!(format!("{:.2}", $b.im)));
		)
	}

	#[test]
	fn forward() {
		let mut output = vec![Complex::new(0.0, 0.0); 4];
		super::forward(Stride::new(&[1.0, 1.0, 0.0, 0.0]), output.as_stride_mut());

		assert_approx_eq!(output[0], Complex::new(2.00,  0.00));
		assert_approx_eq!(output[1], Complex::new(1.00, -1.00));
		assert_approx_eq!(output[2], Complex::new(0.00,  0.00));
		assert_approx_eq!(output[3], Complex::new(1.00,  1.00));
	}

	#[test]
	fn inverse() {
		let mut output = vec![Complex::new(0.0, 0.0); 4];
		super::inverse(Stride::new(&[1.0, 1.0, 0.0, 0.0]), output.as_stride_mut());

		for output in output.iter_mut() {
			ComplexMut::unscale(output, 4.0);
		}

		assert_approx_eq!(output[0], Complex::new(0.50,  0.00));
		assert_approx_eq!(output[1], Complex::new(0.25,  0.25));
		assert_approx_eq!(output[2], Complex::new(0.00,  0.00));
		assert_approx_eq!(output[3], Complex::new(0.25, -0.25));
	}
}