//@code
proc FFT() {local n, x
        if ($1 == 1) { // forward
                $o3.fft($o2, 1)
                n = $o3.size()
                $o3.div(n/2)
                $o3.x[0] /= 2	// makes the spectrum appear discontinuous
                $o3.x[1] /= 2	// but the amplitudes are intuitive

                $o4.copy($o3, 0, 1, -1, 1, 2)   // odd elements
                $o3.copy($o3, 0, 0, -1, 1, 2)   // even elements
                $o3.resize(n/2+1)
                $o4.resize(n/2+1)
                $o3.x[n/2] = $o4.x[0]           //highest cos started in o3.x[1
                $o4.x[0] = $o4.x[n/2] = 0       // weights for sin(0*i)and sin(PI*i)
	}else{ // inverse
                // shuffle o3 and o4 into o2
                n = $o3.size()
                $o2.copy($o3, 0, 0, n-2, 2, 1)
                $o2.x[1] = $o3.x[n-1]
                $o2.copy($o4, 3, 1, n-2, 2, 1)
                $o2.x[0] *= 2
                $o2.x[1] *= 2 
                $o2.fft($o2, -1)
        }
}