function g = czt(x, k, w, a) 
%MATLAB Code Generation Library Function

% Copyright 2018-2023 The MathWorks, Inc.
%#codegen

narginchk(1,4);
% Validate inputs.
% signal X
validateattributes(x,{'double','single'},{'2d'},'czt','X',1)

% transform length M
if nargin < 2
    k = length(x);
end
validateattributes(k,{'numeric'},{'scalar','finite'},'czt','M',2)
% Cast transform length so default complex ratio W is computed in desired
% precision. Input signal X is of type floating point, single or double,
% which is in agreement with data type of W.
k0 = cast(k(1),'like',real(x()));

% ratio between spiral contour points W
if nargin < 3
  w = exp(-1i.*2.*pi./k0); 
end
validateattributes(w,{'double','single'},{'scalar'},'czt','W',3)

% spiral contour initial point A
if nargin < 4
    a = cast(1,'like',coder.internal.scalarEg(x,w));
end
validateattributes(a,{'double','single'},{'scalar'},'czt','A',4)

% Cast inputs to enforce single-precision computation and enforce
% scalarness for scalar inputs W and A.
odt = coder.internal.scalarEg(x,w,a); % output data type
odtr = real(odt);
x1 = cast(x,'like',odt);
k1 = cast(k0,'like',odtr);
w1 = cast(w(1),'like',odt);
a1 = cast(a(1),'like',odt);

% Check whether the input arguments are constant
allConst = true;

allConst = allConst && coder.internal.isConst(x1);
allConst = allConst && coder.internal.isConst(k1);
allConst = allConst && coder.internal.isConst(w1);
allConst = allConst && coder.internal.isConst(a1);

if allConst && coder.internal.isCompiled
    
    % codegeneration for constant input argument
    
    g = coder.const(feval('czt',x1,k1,w1,a1));
    return;
else
    
    % codegeneration for variable input
    % Convert to a column vector if x is a row vector.
    xIsRow =  coder.internal.isConst(isrow(x)) && isrow(x1) && ~isscalar(x1);
    ONE = coder.internal.indexInt(1);
    if xIsRow
        x2 = reshape(x1,[],ONE);
    else
        x2 = x1;
    end
    
    [m, n] = size(x2);

    %------- Length for power-of-two fft.

    nfft = 2^nextpow2(m+k1-1);

    %------- Premultiply data.

    maxDim = coder.internal.indexInt((max(k1-1,m-1) - (-m+1))+1);
    ww = coder.nullcopy(zeros(maxDim,1,'like',odt));
    initVal = cast((-m),'like',odtr)+1;
    for i= 1:maxDim    
        kk = (initVal^2)/2;
        ww(i) = w1^kk;
        initVal = initVal + 1;
    end

    aa= coder.nullcopy(complex(zeros(1,m,'like',odtr),1));
    aa(1) = ww(m);
    for nn = 1:m-1
        aa(nn+1) = (a^-nn)* ww(m+nn);
    end

    y = coder.nullcopy(complex(zeros(m,n,'like',odtr),1));
    for i=1:n
        for j =1:m
            y(j,i) = x2(j,i) .* aa(j);
        end
    end

    %------- Fast convolution via FFT.

    fy = fft(  y, nfft );

    fv = fft( 1 ./ ww(1:(k1-1+m)), nfft );   % <----- Chirp filter.

    for i=1:n
        for j =1:nfft
            fy(j,i) = fy(j,i) .* fv(j);
        end
    end
    
    g1  = ifft( fy );

    %------- Final multiply.
    g2 = coder.nullcopy(complex(zeros(k1,n,'like',odtr),1));
    for i=1:n
        idx = ONE;
         for j=m:m+k1-1
            g2(idx,i) = g1(j,i) .* ww(j);
            idx = idx + 1;
         end
    end

    if xIsRow
        g = reshape(g2,ONE,[]);
    else
        g = g2;
    end
end

end
