function varargout = cpsd(x,y,varargin)
%CPSD Cross power spectral density estimate via Welch's method on GPU.
%
%   Pxy = CPSD(X,Y)
%   Pxy = CPSD(X,Y,WINDOW)
%   Pxy = CPSD(X,Y,WINDOW,NOVERLAP)
%   [Pxy,W] = CPSD(X,Y,WINDOW,NOVERLAP,NFFT)
%   [Pxy,W] = CPSD(X,Y,WINDOW,NOVERLAP,W)
%   [Pxy,F] = CPSD(X,Y,WINDOW,NOVERLAP,NFFT,Fs)
%   [Pxy,F] = CPSD(X,Y,WINDOW,NOVERLAP,F,Fs)
%   [...] = CPSD(...,'mimo')
%   [...] = CPSD(...,FREQRANGE)
%   CPSD(...)
%
%   See also CPSD, GPUARRAY.

%   Copyright 2020-2023 The MathWorks, Inc.

narginchk(2,8);
nargoutchk(0,2);

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

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

onGPU = false;
% Both input signals, window (if it is a nonempty numeric vector), and
% NFFT (if it is a numeric vector) determine whether the calculation is
% done on the GPU and the output datatype
if isgpuarray(x) || isgpuarray(y) || windowOrNFFTOnGPU
    onGPU = true;
    % Ensure that the input signal is on the GPU
    x = gpuArray(x);
    y = gpuArray(y);
    % 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 except for NFFT
    [varargin{extraArgInd}] = gather(varargin{extraArgInd});
else
    % Gather all extra arguments
    [varargin{:}] = gather(varargin{:});
end

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

funcName = 'cpsd';
% Possible outputs are:
%       Plot
%       Pxy
%       Pxy, freq
[varargout{1:nargout}] = welch({x,y},funcName,varargin{:});
end