function varargout = periodogram(x,varargin)
%PERIODOGRAM  Power Spectral Density (PSD) estimate via periodogram method on the GPU.
%
%   Pxx = PERIODOGRAM(X)
%   Pxx = PERIODOGRAM(X,WINDOW)
%   Pxx = PERIODOGRAM(X,WINDOW,...,SPECTRUMTYPE)
%   [Pxx,W] = PERIODOGRAM(X,WINDOW,NFFT)
%   [Pxx,W] = PERIODOGRAM(X,WINDOW,W)
%   [Pxx,F] = PERIODOGRAM(X,WINDOW,NFFT,Fs)
%   [Pxx,F] = PERIODOGRAM(X,WINDOW,F,Fs)
%   [...] = PERIODOGRAM(X,WINDOW,NFFT,...,FREQRANGE)
%   [Pxx,F,Pxxc] = PERIODOGRAM(...,'ConfidenceLevel',P)
%   [RPxx,F] = PERIODOGRAM(X,WINDOW,...,'reassigned')
%   [RPxx,F,Pxx,Fc] = PERIODOGRAM(...,'reassigned')
%   PERIODOGRAM(...)
%
%   See also GPUARRAY, PERIODOGRAM.

%   Copyright 2020-2023 The MathWorks, Inc.

narginchk(1,10);

% Validate input x (validate x early to error out early)
validateattributes(x,{'single','double'},...
    {'nonempty','nonsparse','nonnan','finite'},'periodogram','X');

% Check if WINDOW or NFFT is on GPU
% WINDOW and NFFT are on GPU when each is a vector and a gpuArray
isNFFTSecondArg = true;
[windowOrNFFTOnGPU,windowIsVector,nfftIsVector,nfftArgInd,extraArgInd] = ...
    signal.internal.gpu.utilities.isWindowNFFTOnGPU(isNFFTSecondArg,varargin{:});

onGPU = false;
% Input signal, window (if it is a non-empty numeric vector), and NFFT
% (if it is a non-empty numeric vector) determine whether the
% calculation is done on the GPU and the output datatype.
if isgpuarray(x) || windowOrNFFTOnGPU
    onGPU = true;
    % Ensure that the input signal is on the GPU.
    x = gpuArray(x);
    % Check that a window vector has been provided. If so, ensure that
    % it is on the GPU.
    if windowIsVector
        varargin{1} = gpuArray(varargin{1});
    end
    % If a vector NFFT is provided, make sure it is on GPU
    if nfftIsVector
        varargin{nfftArgInd} = gpuArray(varargin{nfftArgInd});
    end
    % Gather all extra arguments if any 
    [varargin{extraArgInd}] = gather(varargin{extraArgInd});
else
    % Gather all extra arguments.
    [varargin{:}] = gather(varargin{:});
end

if nargin > 1
    [varargin{:}] = convertStringsToChars(varargin{:});
end

% Dispatch in-memory periodogram if input signal is not gpuArray or input
% window vector is not gpuArray.
if ~onGPU
    varargout = cell(1,nargout);
    [varargout{:}] = periodogram(x,varargin{:});
    return;
end

% look for psd, power, and ms window compensation flags
[esttype, args]= signal.internal.psdesttype({'psd','power','ms'},'psd',varargin);
if strcmpi(esttype,'ms')
    error(message("signal:periodogram:InvalidMsLegacyOption"));
end

if isvector(x)
    N = length(x); % Record the length of the data
else
    N = size(x,1);
end

% extract window argument
if ~isempty(args) && ~ischar(args{1})
    win = args{1};
    validateattributes(win,{'double','single'},{'nonsparse'},'periodogram','Window')
    args1 = args(2:end);
else
    win = cast([],'like',real(x()));
    args1 = args;
end

% scan for options
xIsReal = isreal(x);
options = periodogram_options(xIsReal,N,args1{:});

if ~options.reassign && nargout>3
    error(message('MATLAB:nargoutchk:tooManyOutputs'));
end

Fs    = options.Fs;
nfft  = options.nfft;

% Generate a default window if needed
%   default window is rectangular of same size as input.
%   if 'reassigned' is specified, then a Kaiser window is used instead
winName = 'User Defined';
winParam = '';
if isempty(win)
    if options.reassign
        win = kaiser(N,38);
    else
        win = rectwin(N);
        winName = 'Rectangular';
        winParam = N;
    end
    % Ensure that the default window is also on the GPU.
    win = gpuArray(win);
end

if (numel(x)<2 || numel(size(x))>2)
    error(message('signal:computeperiodogram:NDMatrixUnsupported'));
end

% Cast to enforce precision rules. Both x and window are gpuArray.
if superiorfloat(x, win, nfft) == "single"
    x = single(x);
    win = single(win);
    nfft = single(nfft);
    Fs = single(Fs);
end

% Ensure frequency vector is linearly spaced when performing reassignment
checkEquallySpacedFrequencyVectorForReassignment(nfft,options.reassign)

% Compute the PS using periodogram over the whole Nyquist range.
[Sxx,w,RSxx,wc] = computeperiodogram(x,win,nfft,esttype,Fs,options);

% compute reassigned spectrum if needed.
if options.reassign
    RPxx = computepsd(RSxx,w,options.range,nfft,Fs,esttype);
    % Compute the one-sided corrected frequency vector for each component
    % of input signal
    wc = makeOnesidedWc(options.range,wc,nfft);
    % de-alias
    wc = dealiasWc(options.centerdc,Fs,nfft,wc);
else
    RPxx = cast([],'like',RSxx);
    wc = zeros([0,0],'like',w);
end

% Compute the 1-sided or 2-sided PSD [Power/freq] or mean-square [Power].
% Also, compute the corresponding freq vector & freq units.
[Pxx,w,units] = computepsd(Sxx,w,options.range,nfft,Fs,esttype);

% compute confidence intervals if needed.
if ~isnan(options.conflevel)
    % Cast confidence level to the correct precision
    % This casting is done to make sure that computation of confidence
    % interval is performed in the correct precision. Since data type of Fs
    % also indicates whether computation is performed in single precision
    % or double precision, confidence level is cast to data type of Fs.
    cl = cast(options.conflevel,'like',Fs);
    Pxxc = signal.internal.spectral.confInterval(cl, Pxx, xIsReal, w, Fs, 1);
elseif nargout>2 && ~options.reassign
    cl = cast(0.95,'like',Fs);
    Pxxc = signal.internal.spectral.confInterval(cl, Pxx, xIsReal, w, Fs, 1);
else
    Pxxc = cast([],'like',Pxx);
end

if nargout==0 % Plot when no output arguments are specified
    [Pxx,w,Pxxc,RPxx] = gather(Pxx,w,Pxxc,RPxx);
    signal.internal.spectral.plotPeriodogram(Pxx,w,Pxxc,RPxx,options,esttype,units,winName,winParam);
else
    if options.centerdc
        if options.reassign
            RPxx = signal.internal.spectral.psdcenterdc(RPxx, w, [], options);
            wc = centerWc(options,wc);
        end
        [Pxx, w, Pxxc] = signal.internal.spectral.psdcenterdc(Pxx, w, Pxxc, options);
    end
    
    % assign outputs in correct order
    if options.reassign
        % first and third output are reassigned and unnormalized power
        Px = RPxx;
        Px0 = Pxx;
    else
        % first and third output are normalized power and confidence
        % intervals
        Px = Pxx;
        Px0 = Pxxc;
    end
    
    % If the input is a vector and a row frequency vector was specified,
    % return output as a row vector for backwards compatibility.
    if size(nfft,2)>1 && isvector(x)
        Px = real(Px.');
        w = w.';
        if options.reassign
            Px0 = Px0.';
            wc = wc.';
        end
    else
        Px = real(Px);
    end
    
    % Cast to enforce precision rules
    % Only case if output is requested, otherwise plot using double
    % precision frequency vector.
    w = cast(w,'like',Px);
    if options.reassign
        wc = cast(wc,'like',Px);
    else
        wc = zeros(0,0,'like',Px);
    end
    
    outputs = {Px,w,Px0,wc};
    varargout = outputs(1:nargout);
end
end

%------------------------------------------------------------------------------
function options = periodogram_options(isreal_x,N,varargin)
%PERIODOGRAM_OPTIONS   Parse the optional inputs to the PERIODOGRAM function.
%   PERIODOGRAM_OPTIONS returns a structure, OPTIONS, with following fields:
%
%   options.nfft         - number of freq. points at which the psd is estimated
%   options.Fs           - sampling freq. if any
%   options.range        - 'onesided' or 'twosided' psd
%   options.centerdc     - true if 'centered' specified

% Generate defaults
nfftTmp = [];
nfftIsVector = false;
for i=1:length(varargin)
    if isnumeric(varargin{i})
        if (i > 1) && matlab.internal.datatypes.isScalarText(varargin{i-1}) && strncmpi(varargin{i-1},'ConfidenceLevel',strlength(varargin{i-1}))
            % If the current input argument is numeric but the previous input
            % argument is equal to 'ConfidenceLevel', then this numeric
            % argument is the value argument of the name-value pair
            % ConfidenceLevel=CL. Hence, that numeric value is not a value of
            % NFFT, so continue to the next iteration.
            validateattributes(varargin{i},{'double','single'},...
                {'scalar','real','nonempty','nonsparse','positive','<',1},'periodogram','ConfidenceLevel')
            continue;
        else
            nfftTmp = varargin{i};
            % Validate NFFT
            if numel(nfftTmp) > 1
                % Frequency vector
                validateattributes(nfftTmp,{'double','single'},{'real','vector','nonsparse'},'periodogram','F')
                nfftIsVector = true;
            else
                % number of FFT points (NFFT)
                validateattributes(nfftTmp,{'numeric'},{'integer','positive','nonsparse'},'periodogram','NFFT')
            end
            break;
        end
    end
end

if nfftIsVector
    options1.nfft = nfftTmp;
else
    options1.nfft = max(256, 2^nextpow2(N));
end

%assign invalid default value to infer size as 1*1 for codegen
options1.Fs = nan; % Work in rad/sample
options1.centerdc = false;
options1.range = 'twosided';  %Two Sided Range
options1.conflevel = nan; % Default Invalid value

% extract and remove any 'reassigned' option from the argument list
[reassign,args1] = getReassignmentOption(varargin{:});
options1.reassign = reassign;

% Set frequency range for real signals if frequency is not a vector
if isreal_x && ~nfftIsVector
    options1.range = 'onesided'; %One Sided Range
end

% Scan legacy freqrange options and error
if any(strcmp(args1, 'whole'))
    error(message('signal:periodogram:InvalidFreqrangeLegacyOption'));
elseif any(strcmp(args1, 'half'))
    error(message('signal:periodogram:InvalidFreqrangeLegacyOption'));
end

[options,msg,msgobj] = psdoptions(isreal_x,options1,args1{:});

if ~isempty(msg)
  error(msgobj);
end

% Issue a warning and set the range option to 'twosided' if frequency
% vector is specified along with options 'onesided' or 'centered'
if nfftIsVector
    % Frequency vector was specified with 'onesided' or 'centered', return
    % and plot two-sided PSD
    if strcmpi(options.range,'onesided')
        coder.internal.warning('signal:periodogram:InconsistentRangeOption','onesided');
    elseif options.centerdc
        coder.internal.warning('signal:periodogram:InconsistentRangeOption','centered');
        options.centerdc = false;
    end
    options.range = 'twosided';
end

end

%--------------------------------------------------------------------------
function checkEquallySpacedFrequencyVectorForReassignment(nfft,reassign)
% ensure frequency vector is linearly spaced when performing reassignment
if reassign && ~isscalar(nfft)
    f = reshape(nfft,[],1);
    
    % see if we can get a uniform spacing of the freq vector
    [~, ~, ~, maxerr] = signal.internal.spectral.getUniformApprox(f);
 
    if maxerr >= 3*eps(underlyingType(f))
        error(message("signal:periodogram:ReassignFreqMustBeUniform"));
    end
end

end

%--------------------------------------------------------------------------
function [reassign, args]  = getReassignmentOption(varargin)
% search for 'reassigned' flag and return true if present.
% error if 
reassign = false;
args = varargin;
matchLoc = find(strncmpi('reassign',args,8));
if ~isempty(matchLoc)
    reassign = true;
    if any(strcmpi('confidencelevel',args))
        %required to be a runtime error for function signatures : func(x,'conflevel',c)
        error(message('signal:periodogram:ConflictingOptions','reassigned', 'ConfidenceLevel'));
    end
end

args(matchLoc) = [];

end

%--------------------------------------------------------------------------
function w = makeOnesidedWc(range,w,nfft)

if strcmp(range,'onesided') && isscalar(nfft)
   if mod(nfft,2) ~= 0
      select = matlab.internal.ColonDescriptor(1,1,(nfft+1)/2);  % ODD
   else
      select = matlab.internal.ColonDescriptor(1,1,nfft/2+1);    % EVEN
   end
   w = w(select,:);
end
end

%--------------------------------------------------------------------------
function Wc = centerWc(options,Wc)

nFreq = size(Wc,1);
if options.centerdc
  nfft_even = signalwavelet.internal.iseven(options.nfft);
  if strcmp(options.range,'onesided')
    if nfft_even
      Wc = [-Wc(matlab.internal.ColonDescriptor(nFreq-1,-1,2), :) ; ...
          Wc(matlab.internal.ColonDescriptor(1,1,nFreq), :)];
    else
      Wc = [-Wc(matlab.internal.ColonDescriptor(nFreq,-1,2), :) ; ...
          Wc(matlab.internal.ColonDescriptor(1,1,nFreq), :)];
    end
  else % two-sided
    if nfft_even
      Wc = [Wc(matlab.internal.ColonDescriptor(nFreq/2+2,1,nFreq), :) ; ...
          Wc(matlab.internal.ColonDescriptor(1,1,nFreq/2+1), :)];
    else
      Wc = [Wc(matlab.internal.ColonDescriptor((nFreq+1)/2+1,1,nFreq), :) ; ...
          Wc(matlab.internal.ColonDescriptor(1,1,(nFreq+1)/2), :)];
    end
  end
end
end

%--------------------------------------------------------------------------
function Wc = dealiasWc(centerdc,Fs,nfft,Wc)
if isnan(Fs)
    % normalize to 2*pi when default specified.
    Fs1 = cast(2*pi,'like',Fs);
else
    Fs1 = Fs;
end

hasNegativeFreqs  = any(nfft < 0);

if centerdc || hasNegativeFreqs
   % map to [-Fs/2,Fs/2) when using negative frequencies
   Wc = mod(Wc+Fs1/2,Fs1)-Fs1/2;
else
   % map to [0,Fs) when using positive frequencies
   Wc = mod(Wc,Fs1);
end
end