function [amp, specAmp, freq, harmSum, totalNoise, maxSpur] = HarmonicDistortion(xdata, ydata, rbw, isPowerSpec, isDecibelUnit, isTwoSided, numHarmonics) 
% Compute harmonic distortion metrics
%   xdata         - frequency input vector
%   ydata         - power spectrum / PSD input vector
%   isDecibelUnit - true if ydata in units of dB (Watts assumed otherwise)
%   isTwoSided    - true if x/y data represents a two-sided (complex) spectrum
%   numHarmonics  - number of harmonics to compute (for SNR/THD)

% This file is for internal use only.

%   Copyright 2013 The MathWorks, Inc.
%#codegen

% pre-initialize size of freq/amp table
freq = NaN(numHarmonics,1);
amp = NaN(numHarmonics,1);

% convert to one-sided power density
[xdata, ydata] = computeOneSidedPSD(xdata, ydata, rbw, isPowerSpec, isDecibelUnit);

if isempty(ydata)
  specAmp = amp;
  harmSum = NaN;
  totalNoise = NaN;
  maxSpur = NaN;
  return
end

df = mean(diff(xdata));

numCleared = 0;

% DC component may be prominent and leak into low frequency bins
% remove DC component monotonically decreasing from both directions
idxDC = find(abs(xdata)<=df/2,1,'first');
if ~isempty(idxDC)
  if ~isTwoSided
    ydata(idxDC) = 2*ydata(idxDC);
  end
  [ydata, numDC] = clearLocalPeak(ydata, idxDC);
  numCleared = numCleared + numDC;
end

% obtain a rough estimate of the frequency of the maximum component
% we will use this our first harmonic (fundamental) frequency
[f, freq(1), amp(1), iMaxAmp] = computeMaxFrequency(xdata, ydata, rbw, df);

% remove the content of the peak so that it does not interfere with
% non-harmonic (spurious and noise) computation
[ydata, numFund] = clearLocalPeak(ydata, iMaxAmp);
numCleared = numCleared + numFund;
harmSum = 0;
for i=2:numHarmonics
  % find closest bin to i * fundamental freq
  iBin = find(abs(xdata-i*f)<df/2,1,'first');   
  % obtain local maximum value in neighborhood of bin        
  if ~isempty(iBin)
    iLeftBin = max(1,iBin-1);
    iRightBin = min(iBin+1,numel(ydata));
    [ampHarm, idxHarm] = max(ydata(iLeftBin(1):iRightBin(1)));
    freq(i) = xdata(iLeftBin+idxHarm-1);
    amp(i) = rbw*ampHarm;
    harmSum = harmSum + amp(i);
    % remove tails of peak
    [ydata, numHarm] = clearLocalPeak(ydata, iLeftBin+idxHarm-1);
    numCleared = numCleared + numHarm;
  end
end

remainingDistortion = sum(ydata) * df;

% largest spur excluding harmonic content
maxSpur = max(ydata) * rbw;

% largest spur including harmonic content
maxSpur = max([maxSpur; amp(2:end)]);
estimatedNoiseDensity = median(ydata(ydata>0));

% if more than 70% of signal is due to window leakage
% do not return a noise estimate
if numCleared > 0.7 * numel(xdata)
  totalNoise = NaN;
else
  totalNoise = remainingDistortion + numCleared * df * estimatedNoiseDensity;
end

specAmp = computeSpectrumFromPSD(amp,rbw,isPowerSpec,isDecibelUnit);