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