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

14 visualizzazioni (ultimi 30 giorni)
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()?

Risposta accettata

David Goodmanson
David Goodmanson il 21 Set 2019
Modificato: David Goodmanson il 21 Set 2019
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))
.
  1 Commento
Peter Maxwell
Peter Maxwell il 21 Set 2019
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.

Accedi per commentare.

Più risposte (1)

Ajay Pattassery
Ajay Pattassery il 20 Set 2019
fftshift will shift the zero frequency components to the center of the spectrum.
  1 Commento
Peter Maxwell
Peter Maxwell il 20 Set 2019
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.

Accedi per commentare.

Categorie

Scopri di più su Language Fundamentals in Help Center e File Exchange

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by