If you consider ideal conditions in the measurements (anechoic conditions - no reflections, ideal omnidirectionality of the microphones), then you could use the recorded magnitude to represent the polar plots of your source's directivity. Most probably you would like to do it in 2D polar plots for various azimuths and elevations. You can also choose to plot either vertical or horizontal directivity patterns in this way. Alternatively, you could choose to create a "balloon plot" like the ones given by the Common Loudspeaker Format (CLF) files.
Now, let's briefly describe the process for a 2D horizontal polar plot. In the easy case where you have $N$ microphones lying on the same horizontal plane, you just plot the recorded magnitude of each of them and find a way to fill in the gaps. This means that in the simplest case where they are spaced in equal angle steps around the source you could just interpolate every angle between them (by the way, if I am not mistaken, this is what happens in the CLF files, as their angular step is around 5 degrees). Keep in mind that this is NOT the most accurate/best way to do it.
Now, if you happen to either have a non-constant angular step you would have to find an appropriate formula to calculate the magnitude in between the microphone positions. In addition to that, if you happen to plot directivity patterns of horizontal planes that you either have some microphones on for some angles, but not for all of them (i.e. 3 mics for some horizontal angles lie on the plane but another $N - 3$ mics don't lie on the plane) you would also have to interpolate between adjacent microphone positions to get a magnitude value for that angle. A somewhat simple way to do it would be to use Vector Base Amplitude Panning (VBAP) (although this algorithm is used for "panning" it could very well serve as the basis for other calculations) to calculate the linear contribution of adjacent microphones to each location. Nevertheless, you would eventually have to find a way to calculate the magnitude for a $p = [x, y, z]^{T}$ point in 3D space in order to plot the magnitude on a $P = [x, y, c]^{T}$ or $P = [c, y, z]^{T}$ plane.
Now, regarding the frequency representation of your polar plots. If you consider the result of the FFT as a filter bank of linearly spaced equal bandwidth filters, then you can treat the data as a filter-banked data set and do the plotting for each frequency bin separately. This would result in most cases with a humongous amount of data and plots which wouldn't really be any useful (especially for high frequencies the human auditory system shows decreased frequency discrimination. For more info see critical bands). If you would like to have data in a "more meaningful human-related" way (this is somewhat a heavy statement here, take it a bit light though) you could use filterbanks in either 1/3rd-octave spaced bands (which is quite common), 1/6th-octave spaced bands (less common but a bit closer to the human auditory system pitch JND and critical bands bandwidth) or even use something like the bark scale or the mel scale to space your filters. Of course, the same applies to balloon plots.
Regarding your question about the way to "average" frequency bins to create "fractional-octave" plots, keep in mind that simple averaging, or even power/energy summation of the bins doesn't necessarily correspond to the intended "fractional octave" band of interest. For more information, you could have a look at "Generalized Fractional-Octave Smoothing of Audio and Acoustic Responses" by Hatziantoniou and Mourjopoulos.
