Finding the right measure to compare sound signals in the frequency domain
[EDIT: a related discussion on (power) spectrum centroids can be found in How to calculate the mean/center frequency of the spectrum?]

Let us assume you have a relatively concentrated peak, with noise. I assume your maximum frequency is ok. If too high, you can trim the spectrum to lower frequency, which could help. 

The code below gives a simple model, estimates a peak location, and a gross measure of dispersion (or uncertainty). It is based on a weighted average or center of mass. If the peak is narrow, the center of mass (and the median) will be drifted away toward the center of the frequency axis ($f_i$'s). Because even if low, the noise spreads. This is an illustration of Archimedes' lever principle:


  Give me a place to stand on, and I can move the earth.




To achieve your goal, you can give more weight to your spectrum of amplitude $a_i$, and compute an estimated peak location $\bar{f}$:
$$\bar{f} = \frac{\sum_i w_i f_i}{\sum_i w_i}\,. $$

If you choose $w_i = a_i$, you get the standard average. This is not enough. You can increase the peak weight with a power: $w_i = a_i^p$. This is the aim of weightPower = 4 in the code. The location is given by the red o, and the x are the left and right bounds derived from a weighted standard deviation. The higher the $p$, the sharper the estimation, with a risk of major errors if the data differs from the model.

Please report your tests, to adjust the answer. And the result follows. With $p=1$, you get $360$, with $p=2$, you get $243$, with  $p=4$, $\bar{f}=201$, quite close to $200$ in my model. If you need much higher robustness, you can go for weighed medians instead of weighed means.



freqAxis = linspace(0,1000,1000)';
freqPeak = 200;
freqPeakAmplitude = 1;
freqGGDExponent = 2;
freqGGDWidth = 50;
freqNoise = 0.2;

dataSpectrum = exp(-((freqAxis-freqPeak)/freqGGDWidth).^freqGGDExponent)+freqNoise*rand(size(freqAxis));

weightPower = 4; 
freqEstLocation = sum((dataSpectrum.^weightPower).*freqAxis)/sum(dataSpectrum.^weightPower);
freqEstDispersion = sqrt(sum((dataSpectrum.^weightPower).*(freqAxis-freqEstLocation).^2)/sum(dataSpectrum.^weightPower));

freqEstAmplitude = interp1(freqAxis,dataSpectrum,freqEstLocation);
clf;hold on
plot(freqAxis,dataSpectrum);
plot(freqEstLocation,freqEstAmplitude,'or');
plot([freqEstLocation-freqEstDispersion freqEstLocation+freqEstDispersion],freqEstAmplitude,'xr');
axis tight;grid on