OTFS Implementation (Twisted Convolution ) is that a vaild implementation for the twisted Convolution?

function Y = TwistedConv(H, X)
% TwistedConv  Discrete "twisted" convolution in the Delay-Doppler domain
%              (incl. handling of quasi-periodicity)
%
%   Y = TwistedConv(H, X)
%
% Inputs:
%   H : MxN complex array (Channel h_dd)
%   X : MxN complex array (Input x_dd)
%
% Output:
%   Y : MxN complex array (Output y_dd)
%
% Formula (implemented):
%   Y[k',l'] = sum_{k=0}^{M-1} sum_{l=0}^{N-1}
%              H[k,l] * X[(k'-k)_M, (l'-l)_N] * exp(j*2*pi*(k'-k)*l/(M*N))
%
% Note:
%   MATLAB indices are 1-based -> +1 and mod(...)+1 used in the code.

[M, N] = size(H);
assert(isequal(size(X), [M, N]), 'H and X must have the same size MxN.');

Y = zeros(M, N);

for kp = 0:M-1
for lp = 0:N-1
    acc = 0;
    for k = 0:M-1
        for l = 0:N-1
            % raw differences
            dk = kp - k;
            dl = lp - l;
            % We bring dk, dl into the range [0, M-1] and [0, N-1] respectively
            % Note: the formula uses the modulo index in X, which is handled
            % by the MATLAB indexing X(mod(dk,M)+1, mod(dl,N)+1).
            % However, the quasi-periodicity phase factors below *depend* on
            % the un-wrapped differences (or equivalently, on the indices
            % outside the [0, M-1]/[0, N-1] range *before* the modulo
            % operation implied by the 1-based indexing into the MxN array).
            % We compute the modulo result *only* for the index in X, as per
            % the standard formula's definition of the modulo index.
            % The phase_qp terms below are *not* strictly part of the standard
            % twisted convolution for *periodic* input, but are added here
            % for *quasi-periodic* context (e.g., in Delay-Doppler domain
            % for OQAM/OTFS systems).
            
            % The standard modulo for X index:
            dk_mod = mod(dk, M);
            dl_mod = mod(dl, N);


            % --- Quasi-Periodicity Phases (if crossing period boundary) ---
            % If dk, dl is outside [0, M-1] and [0, N-1] -> Period jump
            % The original code seems to have an error in applying
            % quasi-periodicity correction, as it uses the raw dk/dl
            % but computes phase_k/phase_l only for a single period jump.
            % *For the sake of a correct translation of the given code,*
            % *I translate the logic as it is written,*
            % *even though the quasi-periodicity logic seems incomplete*
            % *(it does not account for multiple period wraps)*.
            
            % The original logic assumes the indices (k'-k) and (l'-l) are
            % used to determine the *quasi-periodicity* phase factor,
            % while the *data* is accessed via the *periodic* index
            % (dk_mod, dl_mod).

            phase_k = 1;
            if dk >= M
                % one period further in Delay
                phase_k = exp(1j * 2*pi * l / N);
            elseif dk < 0 % && dk <= -1 is redundant
                % one period previous in Delay
                phase_k = exp(-1j * 2*pi * l / N);
            end

            phase_l = 1;
            if dl >= N
                % one period further in Doppler
                phase_l = exp(-1j * 2*pi * k / M);
            elseif dl < 0 % && dl <= -1 is redundant
                % one period previous in Doppler
                phase_l = exp(1j * 2*pi * k / M);
            end

            phase_qp = phase_k * phase_l;

            % Corrected "twist" phase factor
            % Note: the formula uses (k'-k) and the code uses dk = kp - k.
            phase_dd = exp(1j * 2*pi * dk * l / (M * N));

            % Summand
            % Use the modulo indices (dk_mod, dl_mod) for accessing X as per the formula X[(k'-k)_M, (l'-l)_N]
            acc = acc + H(k+1, l+1) * X(dk_mod+1, dl_mod+1) * phase_dd * phase_qp;
        end
    end
    Y(kp+1, lp+1) = acc;
end
end
end