function varargout = mscohere(x,y,varargin)
%MSCOHERE Magnitude squared coherence estimate via Welch's method on GPU.
%
%   Cxy = MSCOHERE(X,Y)
%   Cxy = MSCOHERE(X,Y,WINDOW)
%   Cxy = MSCOHERE(X,WINDOW,NOVERLAP)
%   [Cxy,W] = MSCOHERE(X,WINDOW,NOVERLAP,NFFT)
%   [Cxy,W] = MSCOHERE(X,WINDOW,NOVERLAP,W)
%   [Cxy,F] = MSCOHERE(X,WINDOW,NOVERLAP,NFFT,Fs)
%   [Cxy,F] = MSCOHERE(X,WINDOW,NOVERLAP,F,Fs)
%   [...] = MSCOHERE(...,'mimo')
%   [...] = MSCOHERE(...,FREQRANGE)
%   MSCOHERE(...)
%
%   See also MSCOHERE, 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 nonempty 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 mscohere if input signal is not gpuArray or input window
% vector is not gpuArray
if ~onGPU
    varargout = cell(1,nargout);
    [varargout{:}] = mscohere(x,y,varargin{:});
    return;
end

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

if nargout == 0
    title(getString(message('signal:dspdata:dspdata:CoherenceTitle')));
    if any(strcmpi(varargin,'MIMO'))
        ylabel(getString(message('signal:dspdata:dspdata:MultipleCoherence')));  
    else
        ylabel(getString(message('signal:dspdata:dspdata:MagnitudeSquaredCoherence')));  
    end
end
end