function [Pxx,F,RPxx,Fc] = computeperiodogram(x,win,nfft,esttype,Fs,options)
%COMPUTEPERIODOGRAM   Periodogram spectral estimation.
%   This function is used to calculate the Power Spectrum Sxx, and the
%   Cross Power Spectrum Sxy.
%   Input X can be a matrix or a multidimensional array.

%   Copyright 2019-2023 The MathWorks, Inc.

narginchk(5,7);
if nargin<6
  reassign = false;
  range = 'twosided';
else
  reassign = options.reassign;
  range = options.range;
end

% use normalized frequencies when Fs is empty
if isempty(Fs) || isnan(Fs)
    Fs = cast(2*pi,'like',Fs);
end

% Validate inputs and convert row vectors to column vectors.
[x,~,y,is2sig,win] = validateinputs(x,win,nfft);

% Window the data
xw = x.*win;

% Compute the periodogram power spectrum [Power] estimate
% A 1/N factor has been omitted since it cancels
[Xx,F] = computeDFT(xw,nfft,Fs);
realXx = real(Xx([]));
F = cast(F,'like',realXx);
if reassign
    dtwin1 = signal.internal.spectral.dtwin(win,Fs);
    xtw = x.*dtwin1;
    [Xxc,~] = computeDFT(xtw,nfft,Fs);
    Fc = arrayfun(@freqCorrect,Xxc,Xx,F);
else
    Fc = zeros([0,0],'like',F);
end

% Evaluate the window normalization constant.  A 1/N factor has been
% omitted since it will cancel below.
if any(strcmpi(esttype,{'ms','power'}))
  if reassign
    if isscalar(nfft)
      U = nfft*(win'*win);
    else
      U = numel(win)*(win'*win);
    end
  else
    % The window is convolved with every power spectrum peak, therefore
    % compensate for the DC value squared to obtain correct peak heights.
    U = sum(win)^2;
  end
else
    U = win'*win;  % compensates for the power of the window.
end

% If two signals are used, we are being called from welch and are not
% performing reassignment.
if is2sig
  y = y.*win;
  [Yy,F] = computeDFT(y,nfft,Fs);
  % Yy can be a single vector or a matrix
  Pxx = (Xx.*conj(Yy))/U;  % Cross spectrum.
else
  % We do a real cast to make Pxx real in codegen.
  Pxx = real(Xx.*conj(Xx)/U);                % Auto spectrum.
end


% Perform reassignment
if reassign
  RPxx = reassignPeriodogram(Pxx, F, Fc, nfft, range);
else
  RPxx = cast([],'like',Pxx);
  Fc = zeros([0,0],'like',F);
end
end

%--------------------------------------------------------------------------
function [fcorr] = freqCorrect(yc,y,f)
% Apply frequency correction from time derivative window

fcorr = -imag(yc ./ y);
if ~isfinite(fcorr)
    fcorr = zeros("like",f);
end
fcorr = fcorr+f;
end

%--------------------------------------------------------------------------
function [x,Lx,y,is2sig,win] = validateinputs(x,win,~)
% Validate the inputs to computexperiodogram and convert row vectors to
% column vectors for backwards compatiblility with R2014a and prior
% releases

% Set defaults and convert to row vectors to columns.
is2sig = false;
win = reshape(win, [], 1);
Lw = length(win);

% Determine if one or two signal vectors was specified.
if iscell(x)
    if length(x) > 1
        y = x{2};
        validateattributes(y,{'single','double'}, {'nonnan','finite'},'periodogram','x');
        if isvector(y)
            y = reshape(y, [], 1);
        end
        is2sig = true;
    end
    x = x{1};
else
    y = zeros([0,0],'like',x);
end

validateattributes(x,{'single','double'}, {'nonnan','finite'},'periodogram','x');
if isvector(x)
    x = reshape(x, [], 1);
end

Lx = size(x,1);

if is2sig
    Ly  = size(y,1);
    if Lx ~= Ly
        coder.internal.error('signal:computeperiodogram:invalidInputSignalLength');
    end
    if size(x,2)~=1 && size(y,2)~=1 && size(x,2) ~= size(y,2)
        coder.internal.error('signal:computeperiodogram:MismatchedNumberOfChannels');
    end
end

% When we reach this code, we've already done the validation for ND-array
% inputs. We still need to check that the input segments have as many
% samples as the given window length, and that the given input doesn't have
% any of the small dimensions set to 0. This happens when the input is of
% size N-by-0, being N greater than 2 and equal to the window length.
sz = size(x);
Lx = sz(1);
coder.internal.errorIf(Lx ~= Lw, 'signal:computeperiodogram:invalidWindow', 'window');
coder.internal.errorIf(any(sz(2:end) == 0), 'signal:computeperiodogram:NDMatrixUnsupported');
end

% -------------------------------------------------------------------------
function RP = reassignPeriodogram(P, f, fcorr, nfft, range)

% for each column input of Sxx, reassign the power additively
% independently.

nChan = size(P,2);

nf = numel(f);
fmin = head(f,1);
fmax = tail(f,1);

% compute the destination row for each spectral estimate
% allow cyclic frequency reassignment only if we have a full spectrum
fcorrVector = reshape(fcorr, [], 1);
if isscalar(nfft) && strcmp(range,'twosided')
    rowIdx = 1+mod(round((fcorrVector-fmin)*(nf-1)/(fmax-fmin)),nf);
else
    rowIdx = 1+round((fcorrVector-fmin)*(nf-1)/(fmax-fmin));
end

% compute the destination column for each spectral estimate
chanIdx = gpuArray.colon(1, nChan);
colIdx = repmat(chanIdx,nf,1);

% reassign the estimates that fit inside the frequency range
P = reshape(P, [], 1);
idx = find(rowIdx>=1 & rowIdx<=nf);
RP = accumarray([rowIdx(idx) colIdx(idx)], P(idx), [nf nChan]);
end

% -------------------------------------------------------------------------

% LocalWords:  Sxx Sxy NFFT ESTTYPE RSxx Fc Fs Yy computexperiodogram
% LocalWords:  compatiblility nonnan
