Matlab IFFT, wavenumber definitions and scalings

Hello everyone,
I am having some trouble with Matlab's FFT and IFFT implementations. My end goal is to be able to numerically construct a function in spectral space (F(k,z)) and then use IFFT to take the result back to physical space with proper scaling (f(x,z)) when it is known that f(x,z) is a purely real function.
Before charging into my full development, I wanted to test that I understood how to use IFFT properly by first inputting a known Fourier transform (F(k) = 2*alpha/(alpha^2 + k^2)) and making sure that it gave me the correct known physical function (f(x) = exp(-alpha* abs(x)). Note this transform pair follows the pure mathematics scalings of Fourier transforms (see here http://uspas.fnal.gov/materials/11ODU/FourierTransformPairs.pdf).
The IFFT of F(k) does not match f(x) and the energies are completely wrong. So I know I am doing something wrong, but after reading a lot of Matlab Q&A about IFFT, I am still lost. I would certainly appreciate any help!
Thank you!
***************
close all
clear all
%physical space
Nx = 128;
x = linspace(-2, 5, Nx);
dx = x(2)-x(1);
%spectral space, positive values only
k1 = 2*pi/(Nx) *[0:Nx/2]; %in radian/sec
k2 = 1/(dx*Nx) *[0:Nx/2]; %in Hertz
%Functional solution
alpha = 0.3;
f = @(x) exp(-alpha*abs(x));
F = @(k) 2*alpha./(alpha^2 + k.^2);
f = f(x);
F1 = F(k1);
F2 = F(k2);
figure(1)
plot(x, f, 'r');
figure(2);
plot(k1, F1, 'r', k2, F2, 'b--');
F1 = cat(2, F1, conj(F1(end-1:-1:2))); %Possible because Hermitian
f1 = ifft(F1, [],1)*1/(dx*Nx); % Scaling?
F2 = cat(2, F2, conj(F2(end-1:-1:2))); %Possible because Hermitian
f2 = ifft(F2, [],1)*1/(dx*Nx); % Scaling?
figure(1)
hold on
plot(x, f1, 'b--', x, f2, 'b.-')
%energy check
energy_f = sum(f.*conj(f) *dx)
energy_F1 = sum(F1.*conj(F1)*1/(dx*Nx))
energy_F2 = sum(F2.*conj(F2)*1/(dx*Nx))
energy_f1 = sum(f.*conj(f) *dx)
energy_f2 = sum(f.*conj(f) *dx)

1 Comment

P.S. Taking:
x = linspace(-a, a, Nx); %some symmetric range around 0
f1 = ifftshift(ifft(F1, [],1)*1/(dx*Nx));
f2 = ifftshift(ifft(F2, [],1)*1/(dx*Nx));
shows a little more clearly that f1 and f2 are peaked as f is supposed to be, but the comparison emphasizes just how different the peaks actually are.
Thanks!

Sign in to comment.

Answers (0)

Asked:

ACE
on 10 Aug 2013

Community Treasure Hunt

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

Start Hunting!