function [amp, specAmp, freq, upperToi, lowerToi] = ImodDistortion(xdata, ydata, rbw, isPowerSpec, isDecibelUnit, isTwoSided)
% Compute intermodulation 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

% This file is for internal use only.

%   Copyright 2013 The MathWorks, Inc.
%#codegen

%pre-initialze freq and amp
freq = NaN(4,1);
amp = NaN(4,1);
f = zeros(4,1);

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

% leave if no data remains
if isempty(ydata)
  specAmp = amp;
  upperToi = NaN;
  lowerToi = NaN;
  return
end

df = mean(diff(xdata));

% 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 = clearLocalPeak(ydata, idxDC);
end

% search line data for largest frequency 
[f(1), freq(1), amp(1), iMaxAmp] = computeMaxFrequency(xdata, ydata, rbw, df);      
ydata = clearLocalPeak(ydata, iMaxAmp);

% Now find the next largest frequency
[f(2), freq(2), amp(2), iMaxAmp] = computeMaxFrequency(xdata, ydata, rbw, df);      
ydata = clearLocalPeak(ydata, iMaxAmp);

% ensure lower and upper tones are in order
if freq(1) > freq(2)
  amp([1 2]) = amp([2 1]);
  freq([1 2]) = freq([2 1]);
  f([1 2]) = f([2 1]);
end

% find closest bin to lower imod
iBin = find(abs(xdata-2*f(1)+f(2))<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(3) = xdata(iLeftBin+idxHarm-1);
  amp(3) = rbw*ampHarm;
end

% find closest bin to upper imod
iBin = find(abs(xdata-2*f(2)+f(1))<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(4) = xdata(iLeftBin+idxHarm-1);
  amp(4) = rbw*ampHarm;
end

% return spectrum for peak labels
specAmp = computeSpectrumFromPSD(amp,rbw,isPowerSpec,isDecibelUnit);

% always report amp in dBW / dBm
amp = 10*log10(amp);

upperToi = amp(2) + (amp(1)-amp(4))/2;
lowerToi = amp(1) + (amp(2)-amp(3))/2;

