MDCT via FFT

mdct_via_fft.m

function y = mdct_via_fft(x)
    % This was edited together from multiple fns in the original source
    % for this gist, I might well have introduced errors along the way.
    
    % Rows of x correspond to samples
    nrows = size(x,1);
    if mod(nrows,4) ~= 0 || nrows<8
        error('Number of rows must be multiple of 4 and at least 8.');
    end
    nr4 = nrows/4;
    
    % Chop into 4 pieces
    a = x(0*nr4+1:1*nr4,:);
    b = x(1*nr4+1:2*nr4,:);
    c = x(2*nr4+1:3*nr4,:);
    d = x(3*nr4+1:4*nr4,:);
    % Rearranged input vector
    t = [-flipud(c)-d; a-flipud(b)];
    
    % Then compute as DCT-IV of rearranged input vector
    % as per R. Gluth, "Regular FFT-related transform kernels for DCT/DST-based polyphase filter banks" (1991)
    
    % First, interleave into complex values
    u_prime = t(1:2:end,:) + j*flipud(t(2:2:end,:));
    % Calculate twiddles
    twiddle = exp(-j*pi/nrows*((0:nr4-1)' + 1/8));
    % The main calc
    f = twiddle .* fft(twiddle .* u_prime);
    % Build the result
    y = zeros(nr4*2,size(x,2));
    y(1:2:end,:) = real(f);
    y(2:2:end,:) = -flipud(imag(f));
end