MATLAB: Centered FFT & DFT: cannot devise required phase shift vector(s).

centered fftfft phase shiftgdftgeneralized dftshifted fft

Can anyone help with using MATLAB builtin fft() function to do a centered FFT/DFT, see Wikipedia – Discrete Fourier Transform, Generalised / Centered.
The background is that I'm doing FFTs and inverse FFTs between a physical coordinate space and k-space for modelling of dispersive waves. So there are two reasons for requiring the centered FFT:
i. both domains have negative coordinates, centered around the origin; and,
ii. the origin in k-space is often undefined and so I'd like to body-swerve that whole issue by using even-length vectors/arrays, so that the origin itself is not included.
So far, I can do the calculation correctly if I manually construct a DFT matrix and premultiply the vectors (it's ok but less than ideal due to the computational complexity). However, I cannot seem to get this right when using phase shifts with fft().
The shifts I derived are obviously wrong, producing nonsense. I also tried the code here: https://se.mathworks.com/matlabcentral/fileexchange/55929-centered-fourier-transforms The real component is off-by-one and the imaginary component is nonzero when it shouldn't be.
I'm using a Guassian and a characteristic function for test purposes. Attached should be some sample code.
Does anybody know how to do the phase shifts correctly so that it'll work with MATLAB's builtin fft()?

Best Answer

Hi Peter,
here are some phase shifts. I shortened up some of the variable names to make it easier (for me, anyway) to follow what was going on.
The code has an alternate way to calculate the initial dft matrix without for loops [ instead of
exp(2*pi*i/N) ^ (...)   the expression uses    exp((2*pi*i/N)*(...))
since I think it's better. That's just a detail ].
It doesn't really matter in this case but I used 400 points because I think unless it's unavoidable, using 2^n points is almost always a useless and error-prone waste of time. So I get to mention it again here.
The code shows that the dft method and fft method give the same result, and I will leave any plotting to you (the frequency grid scaling may not be correct, but that's a different issue).
N = 400;
width_x = sqrt(N);
dx = width_x / N;
sigma2 = 0.5 * pi/9;
% for centering, toss in some differences from -(N-1)/2 as a check;
% if it works here it works when differences are zero 
a = -(N-1)/2 + .1;
b = -(N-1)/2 - .4;
% Manually construct centered DFT matrix
% with for loops
aWcen_forloop = ones(N,N);               % a_W_centered
Nthroot = exp(-2i*pi/N );
for ii=1:N
    for jj=1:N        
        aWcen_forloop(ii,jj) = Nthroot^( (ii-1+a) * (jj-1+b) );                
    end
end
% all at once
[iiM jjM] = ndgrid(0:N-1,0:N-1);
aWcen = exp((-2*pi*i/N)*(iiM+a).*(jjM+b));
% check, difference is basically 0
max(max(abs(aWcen_forloop - aWcen)))
% Create x and kx vectors
v_x_cen = ( linspace( 0, width_x, N ) - width_x/2 ).';
v_kx_cen = v_x_cen;
% Create initial Gaussian
v_gauss = exp( -(0.5/sigma2) * v_x_cen.^2 );
% Calculate centered transforms of Gaussian
% by dft matrix
v_gauss_dft_fs_cen = (1/sqrt(N)) * aWcen * v_gauss;
% by fft
vpre =  exp(-2*pi*i*(0:N-1)'*a/N);
vpost = exp(-2*pi*i*(0:N-1)'*b/N);
A = exp(-2*pi*i*a*b/N);            % scalar factor
v_gauss_fft_fs_cen = (1/sqrt(N))*A*vpost.*fft(vpre.*v_gauss);
% check, diference is basically 0
max(abs(v_gauss_fft_fs_cen - v_gauss_dft_fs_cen))
.
Related Question