An fft should return a series of complex numbers (could either be rectangular coordinates, or polar: phase and magnitude) for a specific range of frequencies...
I'm still working through the expressions, but I'll bet dollars to donuts that the x and y arrays are the real (x) and imaginary (y) components of the complex numbers that are the result of the transformation.
The absolute value of the sum of the squares of these two components should be the magnitude of the harmonic component at each frequency (conversion to polar).
If the phase is important for your application, keep in mind that the the FFT (as with any phasor) can either be sine referenced or cosine referenced. I beleive sine is the standard, however.
see:
http://www.mathworks.com/help/matlab/ref/fft.html
http://mathworld.wolfram.com/FastFourierTransform.html
Since the FFT gives a truncated approximation to an infinite series created by a harmonic decomposition of a periodic waveform any periodic waveform can be used to test the functionality of your code.
For an example, a square wave should be easy to replicate, and has very well known harmonic coefficients. The resolution of the data set will determine the number of harmonics that you can calculate (most fft algorithms do best with a data set that has a length equal to a power of two, and is a number of integral wavelengths of the longest frequency that you want to use).
The square wave coefficients should be at odd multiples of the fundamental frequency and have magnitudes that vary inversely with the order of the harmonic.
http://en.wikipedia.org/wiki/Square_wave