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

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()?

Answer 1

[David Goodmanson]

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)) 

.

[Peter Maxwell]

Thanks David, it works perfectly!

Very much appreciated as I was banging my head against a wall the other day with this and just couldn't get it right. It means I can now do the full calculations without requiring a multicore server.

Answer 2

[Ajay Pattassery]

fftshift will shift the zero frequency components to the center of the spectrum.

[Peter Maxwell]

Thanks, Ajay. Yes, I'm aware of fftshit() and ifftshift() but --as far as I know-- no matter whether the number of grid points is even or odd, zero is always included.

So, for even length vectors, one ends up with arrangement akin to

[ -3 -2 -1 0 1 2 ]

rather than the sought

[ -2.5 -1.5 -0.5 0.5 1.5 2.5 ]

See the example .m file I'd included. When using the DFT with the centered shifts, the DFT of the Gaussian distribution has two points that are the same value in the centre whereas when using FFT, there is a single peak value.

Note that both the signal and Fourier domain in the examples I'm working with include negative grid points and are centered around the origin (but, as I said before, I want to explicitly exclude the origin by using even-numbered grid sizes). See also, https://dsp.stackexchange.com/questions/38746/centered-fourier-transform

If I've got this totally wrong then please say but the numerical examples I've tried suggest that I haven't.