Discover MakerZone

MATLAB and Simulink resources for Arduino, LEGO, and Raspberry Pi

Learn more

Discover what MATLAB® can do for your career.

Opportunities for recent engineering grads.

Apply Today

Thread Subject:
Inverse Fourier transforms

Subject: Inverse Fourier transforms

From: Mat Hunt

Date: 17 Jun, 2010 19:34:20

Message: 1 of 32

I need to compute an inverse Fourier transform, I did this previously by writing my own inverse Fourier transform routine which worked pretty well on simple examples but I am currently investigating some rather obscure wave profiles.

I computed these wave profiles two ways and the answers don't agree, so I thought I would delve into the murky depths of the Matlab inverse Fourier transform.

Say for example I have a Fourier transform of the Gaussian function, which is sqrt(pi)*exp(-0.25*k.^2). What commands would I use to get back to my original Gaussian function? I looked on the help and that was of no use, it mentioned the ifft function but not how I should use it.

If I do the following for the Fourier transform k=-100:0.001:100, and the Fourier transform is sqrt(pi)*exp(-0.25*k.^2)

What commands do I use to get back the the original function? I also want to be able to plot it. The function I am using is defined analytically, would it be better to use the symbolic ifourier instead?

Mat

Subject: Inverse Fourier transforms

From: Walter Roberson

Date: 17 Jun, 2010 19:45:30

Message: 2 of 32

Mat Hunt wrote:

> If I do the following for the Fourier transform k=-100:0.001:100, and
> the Fourier transform is sqrt(pi)*exp(-0.25*k.^2)

You appear to be mixing the discrete Fourier transform and the continuous
Fourier transform. The discrete fourier routine is fft() and the inverse
discrete fourier transform is ifft(). I believe the continuous fourier
transform requires the symbolic toolbox (continuous implies infinite if done
numerically) , but perhaps it is supported in some other way as well.

Subject: Inverse Fourier transforms

From: Mat Hunt

Date: 17 Jun, 2010 21:49:04

Message: 3 of 32

I want to compute an inverse Fourier transform and plot the solution. I have the analytical expression of the inverse Fourier transform. I can either supply it as an analytical function (and use ifourier) or give it as a numerical array and use ifft.

Which would be best?

Mat

Subject: Inverse Fourier transforms

From: Walter Roberson

Date: 17 Jun, 2010 21:58:52

Message: 4 of 32

Mat Hunt wrote:
> I want to compute an inverse Fourier transform and plot the solution. I
> have the analytical expression of the inverse Fourier transform. I can
> either supply it as an analytical function (and use ifourier) or give it
> as a numerical array and use ifft.
>
> Which would be best?

I would try ifourier first as it would be more accurate -- using the discrete
fourier is going to give you round-off problems.

If your function is sufficiently complex then there might not be a nice
symbolic inverse transform; in that case you would have to resort to ifft.

Subject: Inverse Fourier transforms

From: Mat Hunt

Date: 17 Jun, 2010 22:29:06

Message: 5 of 32

Walter Roberson <roberson@hushmail.com> wrote in message <hve5un$mj1$1@canopus.cc.umanitoba.ca>...
> Mat Hunt wrote:
> > I want to compute an inverse Fourier transform and plot the solution. I
> > have the analytical expression of the inverse Fourier transform. I can
> > either supply it as an analytical function (and use ifourier) or give it
> > as a numerical array and use ifft.
> >
> > Which would be best?
>
> I would try ifourier first as it would be more accurate -- using the discrete
> fourier is going to give you round-off problems.
>
> If your function is sufficiently complex then there might not be a nice
> symbolic inverse transform; in that case you would have to resort to ifft.

The function in question is:

16.2668*exp(-1/4*k^2)/(144+2*k^4-12*k^2+9*abs(k))

I tried ifourier and I got the result:

16.2668*ifourier(exp(-1/4*k^2)/(144+2*k^4-12*k^2+9*abs(k)),k,x)

So I think I will have to use ifft. I would ideally like to plot the original function but I don't know how to plot it.

Mat

Subject: Inverse Fourier transforms

From: Walter Roberson

Date: 18 Jun, 2010 02:59:39

Message: 6 of 32

Mat Hunt wrote:

> The function in question is:
>
> 16.2668*exp(-1/4*k^2)/(144+2*k^4-12*k^2+9*abs(k))
>
> I tried ifourier and I got the result:
>
> 16.2668*ifourier(exp(-1/4*k^2)/(144+2*k^4-12*k^2+9*abs(k)),k,x)
>
> So I think I will have to use ifft. I would ideally like to plot the
> original function but I don't know how to plot it.

I'm trying to figure out what kind of function could possibly have
abs(k) as part of its fourier transform. It appears that at the very
least the answer involves complex numbers.

Subject: Inverse Fourier transforms

From: Mat Hunt

Date: 18 Jun, 2010 10:04:05

Message: 7 of 32

No, there are no problems with it, no complex numbers. The integrand is an even function and real, as a result, the complex part of the solution is zero.

Mat

Subject: Inverse Fourier transforms

From: Mat Hunt

Date: 18 Jun, 2010 12:58:05

Message: 8 of 32

So, the best way to perform this calculation is numerically. In which case how would I go about this using matlab's inbuilt functions. I wrote a small program myself:

%Start the integration
x=-8:0.01:8;
n=length(x);
y=zeros(n,1);
dk=0.001;
k=-5:dk:5;
l=length(k)-1;
for j=1:n
    z=f(k).*exp(i*k*x(j));
    y(j)=real(trapz(k,z));
end
plot(x,y)
xlabel('x');
ylabel('\eta_{1}(x)');

which seemed to work in some simple cases rather well. Now I want to compare it to matlabs ifft routine.

Any suggestions on what I should do?

Mat

Subject: Inverse Fourier transforms

From: Walter Roberson

Date: 18 Jun, 2010 13:37:34

Message: 9 of 32

Mat Hunt wrote:

> The function in question is:
>
> 16.2668*exp(-1/4*k^2)/(144+2*k^4-12*k^2+9*abs(k))
>
> I tried ifourier and I got the result:
>
> 16.2668*ifourier(exp(-1/4*k^2)/(144+2*k^4-12*k^2+9*abs(k)),k,x)
>
> So I think I will have to use ifft. I would ideally like to plot the
> original function but I don't know how to plot it.

I get something vaguely plausible if I use

k = -100:.1:100;
t = 16.2668*exp(-1/4*k.^2)/(144+2*k.^4-12*k.^2+9*abs(k));
ti = ifft([0 t]);
plot(ti)

Subject: Inverse Fourier transforms

From: Mat Hunt

Date: 18 Jun, 2010 14:18:04

Message: 10 of 32


>
> k = -100:.1:100;
> t = 16.2668*exp(-1/4*k.^2)/(144+2*k.^4-12*k.^2+9*abs(k));
> ti = ifft([0 t]);
> plot(ti)

I used the above code and didn't get anything remotely sensible. The code I use is:

h=1;
E_b=0.1; %0.2;
B=0.2;
rho=1;
F=0.2;
mu=0;
%Start the integration
x=-8:0.01:8;
n=length(x);
y=zeros(n,1);
dk=0.001;
k=-20:dk:20;
l=length(k)-1;
for j=1:n
    z=f_1(k,F,B,E_b,h,rho,mu).*exp(i*k*x(j));
    y(j)=real(trapz(k,z));
end
 
plot(x,y)
xlabel('x');
ylabel('\eta_{1}(x)');

Along with the function:

function eta=f_1(k,F,B,E_b,h,rho,mu)
m=length(k);
for j=1:m
eta(j)=-(sqrt(pi)*exp(-0.25*k(j)^2)/(2*rho*9.8065))./(1-F+(h*k(j))^4/90+0.5*(B-1/3)*(h*k(j))^2-0.5*E_b*h*abs(k(j))+i*mu);
end

If I don't get anything like the above function then the answer is incorrect. The code that you gave me didn't give anything like the above M file.

Mat

Subject: Inverse Fourier transforms

From: Walter Roberson

Date: 19 Jun, 2010 15:53:45

Message: 11 of 32

Mat Hunt wrote:

>> k = -100:.1:100;
>> t = 16.2668*exp(-1/4*k.^2)/(144+2*k.^4-12*k.^2+9*abs(k));
>> ti = ifft([0 t]);
>> plot(ti)

> I used the above code and didn't get anything remotely sensible. The
> code I use is:
>
> h=1;
> E_b=0.1; %0.2;
> B=0.2;
> rho=1;
> F=0.2;
> mu=0;
> %Start the integration
> x=-8:0.01:8;
> n=length(x);
> y=zeros(n,1);
> dk=0.001;
> k=-20:dk:20;
> l=length(k)-1;
> for j=1:n
> z=f_1(k,F,B,E_b,h,rho,mu).*exp(i*k*x(j));
> y(j)=real(trapz(k,z));
> end
>
> plot(x,y)
> xlabel('x');
> ylabel('\eta_{1}(x)');
>
> Along with the function:
>
> function eta=f_1(k,F,B,E_b,h,rho,mu)
> m=length(k);
> for j=1:m
> eta(j)=-(sqrt(pi)*exp(-0.25*k(j)^2)/(2*rho*9.8065))./(1-F+(h*k(j))^4/90+0.5*(B-1/3)*(h*k(j))^2-0.5*E_b*h*abs(k(j))+i*mu);
>
> end

> If I don't get anything like the above function then the answer is
> incorrect. The code that you gave me didn't give anything like the
> above M file.

I missed copying in a "."

t = 16.2668*exp(-1/4*k.^2)/(144+2*k.^4-12*k.^2+9*abs(k));

should have been

t = 16.2668*exp(-1/4*k.^2)./(144+2*k.^4-12*k.^2+9*abs(k));


Anyhow, if you take the fft of your y, you will see that the result is
complex whereas the function you provided earlier is not complex and
thus cannot be the discrete fourier transform of your y.

Is your f_1 intended to be a more detailed version of the function you
gave earlier? If so then note that your f_1 function has an imaginary
component whereas the function you gave earlier does not.

My answer was not incorrect for the question you posted earlier, and I
don't appreciate being expected to read your mind about what you are
really trying to do.

Subject: Inverse Fourier transforms

From: Greg Heath

Date: 19 Jun, 2010 17:10:01

Message: 12 of 32

On Jun 17, 6:29 pm, "Mat Hunt" <hunt_...@hotmail.com> wrote:
> Walter Roberson <rober...@hushmail.com> wrote in message <hve5un$mj...@canopus.cc.umanitoba.ca>...
> > Mat Hunt wrote:
> > > I want to compute an inverse Fourier transform and plot the solution. I
> > > have the analytical expression of the inverse Fourier transform. I can
> > > either supply it as an analytical function (and use ifourier) or give it
> > > as a numerical array and use ifft.
>
> > > Which would be best?
>
> > I would try ifourier first as it would be more accurate -- using the discrete
> > fourier is going to give you round-off problems.
>
> > If your function is sufficiently complex then there might not be a nice
> > symbolic inverse transform; in that case you would have to resort to ifft.
>
> The function in question is:
>
> 16.2668*exp(-1/4*k^2)/(144+2*k^4-12*k^2+9*abs(k))

Are you sure about abs(:) ?. The function is not analytic.

Hope this helps.

Greg

Subject: Inverse Fourier transforms

From: Greg Heath

Date: 19 Jun, 2010 19:12:21

Message: 13 of 32

On Jun 18, 6:04 am, "Mat Hunt" <hunt_...@hotmail.com> wrote:
> No, there are no problems with it, no complex numbers. The integrand is an even function and real, as a result, the complex part of the solution is zero.

No, the function has an imaginary part.

I'll include a demo in my next post.

Greg

Subject: Inverse Fourier transforms

From: Greg Heath

Date: 19 Jun, 2010 19:59:51

Message: 14 of 32

On Jun 18, 9:37 am, Walter Roberson <rober...@hushmail.com> wrote:
> Mat Hunt wrote:
> > The function in question is:
>
> > 16.2668*exp(-1/4*k^2)/(144+2*k^4-12*k^2+9*abs(k))
>
> > I tried ifourier and I got the result:
>
> > 16.2668*ifourier(exp(-1/4*k^2)/(144+2*k^4-12*k^2+9*abs(k)),k,x)
>
> > So I think I will have to use ifft. I would ideally like to plot the
> > original function but I don't know how to plot it.
>
> I get something vaguely plausible if I use
>
> k = -100:.1:100;
> t = 16.2668*exp(-1/4*k.^2)/(144+2*k.^4-12*k.^2+9*abs(k));
> ti = ifft([0 t]);
> plot(ti)

I assume the zero is there to deal with the wild oscillations that
occur in the real and imaginary parts otherwise. I have no idea
why that happens. It also happens when abs(k) is replaced by
k or -k.

The following works well:

close all, clear all, clc, p = 0

kb = -100:.1:100; % bipolar k
dk = kb(2)-kb(1) % 0.1
N = length(kb) % 2001; N is odd
k = dk*(0:N-1); % unipolar k

Yb = 16.2668*exp(-1/4*kb.^2)./(144+2*kb.^4-12*kb.^2+9*abs(kb));
% Yb = 16.2668*exp(-1/4*kb.^2)./(144+2*kb.^4-12*kb.^2+9*kb);
% Yb = 16.2668*exp(-1/4*kb.^2)./(144+2*kb.^4-12*kb.^2-9*kb);
Y = fftshift(Yb);

p = p+1, figure(p) % figure 1
subplot(211)
plot(kb,Yb)
subplot(212)
plot(k,Y)

dx = 1/(N*dk)
x = dx*(0:N-1);
xb = x-dx*(N-1)/2;

y = ifft(Y);
ry = real(y);
iy = imag(y);
ay = abs(y);
py = angle(y);

p = p+1, figure(p) % figure 2
subplot(221)
plot(x,ry)
subplot(222)
plot(x,iy)
subplot(223)
plot(x,ay)
subplot(224)
plot(x,py)

% yb = ifft(Yb);
yb = ifftshift(y);
ryb = real(yb);
iyb = imag(yb);
ayb = abs(yb);
pyb = angle(yb);

p=p+1,figure(p) % figure 3
subplot(221)
plot(xb,ryb)
subplot(222)
plot(xb,iyb)
subplot(223)
plot(xb,ayb)
subplot(224)
plot(xb,pyb)

sortfigs(1,p)

Subject: Inverse Fourier transforms

From: Mat Hunt

Date: 20 Jun, 2010 10:33:06

Message: 15 of 32

Hi Greg,

I guess without knowing the full background to what I am doing. I am working on electrified flow down a channel with a moving pressure distribution. I am working on the weakly nonlinear theory which is a KdV like equation.

I ignored the nonlinear terms and solve the equation via Fourier transforms (the equation contained a Benjamin-Ono term (a Hilbert transform term)) and that is where the |k| term comes in as I am taking the derivative of a Hilbert transform of e^{ikx} which comes out as sgn(k)*k, which is |k|, that is the origin of the |k|.

The complex term is nothing but Rayleigh dispersion, which I have set to zero for the parameters I am using . I solved the equation numerically and I want to compare my answers for the linear case, and currently they don't match at all which is very worrying.

Your solution was close to what I get for my code which is very encouraging.

Mat

Subject: Inverse Fourier transforms

From: Greg Heath

Date: 20 Jun, 2010 19:16:01

Message: 16 of 32

On Jun 20, 6:33 am, "Mat Hunt" <hunt_...@hotmail.com> wrote:
> Hi Greg,
>
> I guess without knowing the full background to what I am doing. I am working on electrified flow down a channel with a moving pressure distribution. I am working on the weakly nonlinear theory which is a KdV like equation.  
>
> I ignored the nonlinear terms and solve the equation via Fourier transforms (the equation contained a Benjamin-Ono term (a Hilbert transform term)) and that is where the |k| term comes in as I am taking the derivative of a Hilbert transform of e^{ikx} which comes out as sgn(k)*k, which is |k|, that is the origin of the |k|.
>
> The complex term is nothing but Rayleigh dispersion, which I have set to zero for the parameters I am using .  I solved the equation numerically and I want to compare my answers for the linear case, and currently they don't match at all which is very worrying.
>
> Your solution was close to what I get for my code which is very encouraging.  

I am intrigued by the rapid oscillations in the real and imaginary
part of ifft(Yb) opposed to the nice and smooth plot of abs(ifft(Yb)).

Hope this helps.

Greg

Subject: Inverse Fourier transforms

From: Mat Hunt

Date: 21 Jun, 2010 08:08:04

Message: 17 of 32

Greg Heath <heath@alumni.brown.edu> wrote in message <4cbcef89-1c16-4d58-8eab-7da8e0e7eff2@r27g2000yqb.googlegroups.com>...
> On Jun 20, 6:33 am, "Mat Hunt" <hunt_...@hotmail.com> wrote:
> > Hi Greg,
> >
> > I guess without knowing the full background to what I am doing. I am working on electrified flow down a channel with a moving pressure distribution. I am working on the weakly nonlinear theory which is a KdV like equation.  
> >
> > I ignored the nonlinear terms and solve the equation via Fourier transforms (the equation contained a Benjamin-Ono term (a Hilbert transform term)) and that is where the |k| term comes in as I am taking the derivative of a Hilbert transform of e^{ikx} which comes out as sgn(k)*k, which is |k|, that is the origin of the |k|.
> >
> > The complex term is nothing but Rayleigh dispersion, which I have set to zero for the parameters I am using .  I solved the equation numerically and I want to compare my answers for the linear case, and currently they don't match at all which is very worrying.
> >
> > Your solution was close to what I get for my code which is very encouraging.  
>
> I am intrigued by the rapid oscillations in the real and imaginary
> part of ifft(Yb) opposed to the nice and smooth plot of abs(ifft(Yb)).
>
> Hope this helps.
>
> Greg

So What do I need to do if I have the Fourier transform F(k), of a function, f as a function of the Fourier transform variable (say k, for example) back to the original variable x, say?

I can send you my other numerical code which I used to compute the integral if you like.

Mat

Subject: Inverse Fourier transforms

From: Greg Heath

Date: 21 Jun, 2010 11:36:43

Message: 18 of 32

On Jun 20, 3:16 pm, Greg Heath <he...@alumni.brown.edu> wrote:
> On Jun 20, 6:33 am, "Mat Hunt" <hunt_...@hotmail.com> wrote:
>
-----SNIP
>
> > Your solution was close to what I get for my code which is very encouraging.
>
> I am intrigued by the rapid oscillations in the real and imaginary
> part of ifft(Yb) opposed to the nice and smooth plot of abs(ifft(Yb)).

% 1. Used trial and error to determine whether to use
% fftshift or ifftshift (different for N odd) for
% shifting from bipolar k to unipolar k:
% Y = ifftshift(Yb)
% and from unipolar x to bipolar x
% y = fftshift(ifft(Y))
% 2. Size is reduced from 2001 to 81 to facilitate
% observing plot details
% 3. Compared yb = fftshift(ifft(Yb)) with y
% a. yb is complex but y is real
% b. They differ by a phase shift that is linear
% in x
% 5. Changing the use of fftshift and ifftshift results in
% both yb and y being complex. The corresponding plots
% for N = 2001 contain relatively wild oscillations (Left
% as an excercise for the reader).

close all, clear all, clc, p = 0

%Original Formulation

kb = -100:.1:100; % bipolar k
N0 = length(kb) % N0 = 2001; N0 is odd
Yb = 16.2668*exp(-1/4*kb.^2)./(144+2*kb.^4-12*kb.^2+9*abs(kb));
minmaxYb = minmax(Yb) % [ 0 0.11296 ]
[Ybmax ibmax] = max(Yb) % [ 0.11296 1001]
[Ymax imax] = max(fftshift(Yb)) % [ 0.11296 2001]
[Ymax imax] = max(ifftshift(Yb)) % [ 0.11296 1 ]

% Convenient Reformulation

kb = -4 : 0.1 : 4; % bipolar k
N = length(kb) % N = 81; N is odd
dk = kb(2)-kb(1) % 0.1
k = dk*(0:N-1); % unipolar k

Yb = 144*exp(-1/4*kb.^2)./(144+2*kb.^4-12*kb.^2+9*abs(kb));
minmaxYb = minmax(Yb) % [ 0.0052749 1 ]
Y = ifftshift(Yb);
[Ymax imax] = max(Y) % [ 1 1 ]

%Similar Transforms

Ybp = 144*exp(-1/4*kb.^2)./(144+2*kb.^4-12*kb.^2+9*kb);
Ybm = 144*exp(-1/4*kb.^2)./(144+2*kb.^4-12*kb.^2-9*kb);

p = p+1, figure(p), hold on % figure 1
plot(kb,Ybp,'r')
plot(kb,Ybm,'g')
plot(kb,Yb)

% Pole placements for analytic inversion
%
%
% Four poles of Yb are at +/-kb0 and +/-conj(kb0)
% No symbolic solution because of abs(k)
% Consider numerical search starting at symbolic
% solution for poles of Ybp and Ybm:

kbp0 = double(solve('144 + 2*kb^4 - 12*kb^2 + 9*kb = 0','kb'))

% kbp0 = 2.3977 - 1.794i
% 2.3977 + 1.794i
% -2.3977 - 1.5099i
% -2.3977 + 1.5099i

kpm0 = -kbp0;

% function denom = denomfunc(kb)
%
% denom = 144+2*kb.^4-12*kb.^2+9*abs(kb)
%
% return

kb0 = fsolve(@denomfunc,kbp0(1),optimset('fsolve'))

% kb0 = 2.4756 - 1.7688i

denomfunc(kb0) % 9.3344e-011 +1.3068e-010i
denomfunc(-kb0) % 9.3344e-011 +1.3068e-010i
denomfunc(conj(kb0)) % 9.3344e-011 -1.3068e-010i
denomfunc(-conj(kb0)) % 9.3344e-011 -1.3068e-010i


% Return to numerical inversion


dx = 1/(N*dk) % 0.12346
x = dx*(0:N-1);
minmax(x) % [ 0 9.8765 ]
xb = x-dx*(N-1)/2;
minmax(xb) % [ -4.9383 4.9383 ]

yb = fftshift(ifft(Yb));
y = fftshift(ifft(Y));
absyb = abs(yb);
absy = abs(y);
maxabs(absy-absyb) % 5.0849e-017
[ymax imax] = max(absy) % [ 0.42304 41]

ryb = real(yb);
ry = real(y);
iyb = imag(yb);
iy = imag(y);

p = p+1, figure(p) % figure 2
subplot(221)
plot(xb,ryb,'b')
title('Real(yb)')
subplot(222)
plot(xb,iyb,'b')
title('Imag(yb)')
subplot(223)
plot(xb,ry,'r')
title('Real(y)')
subplot(224)
plot(xb,iy,'r')
title('Imag(y)')

phaseyb = angle(yb);
phasey = angle(y);

p = p+1, figure(p) % figure 3
subplot(221)
plot(xb,absyb,'b')
title('Mag(yb)')
subplot(222)
plot(xb,phaseyb,'b')
title('Phase(yb)')
subplot(223)
plot(xb,absy,'r')
title('Mag(y)')
subplot(224)
plot(xb,phasey,'r')
title('Phase(y)')

sortfigs(1,p)

Hope this helps.

Greg

Subject: Inverse Fourier transforms

From: Mat Hunt

Date: 21 Jun, 2010 21:42:07

Message: 19 of 32

Hi Greg,

So I need to do an fftshift after the ifft to get the answer? do I also need to do this with the variable? I can evaluate the integral using residue calculus if I ignore the |k| terms and that should provide me with another check on my working.

I had a look at you code and there were some commands which didn't run and that matlab. Are you using an old version?

Mat

Subject: Inverse Fourier transforms

From: Greg Heath

Date: 21 Jun, 2010 23:58:14

Message: 20 of 32

On Jun 21, 5:42 pm, "Mat Hunt" <hunt_...@hotmail.com> wrote:
> Hi Greg,
>
> So I need to do an fftshift after the ifft to get the answer?

Yes, fft and ifft assume the nonnegative intervals x = dx*(0:N-1) and
k = dk*(0:N-1). fftshift and it's inverse ifftshift allow the
translation
to and from the bipolar interval with indexing (0:N-1)-ceil((N-1)/
2).
When N is even, fftshift and ifftshift yield the same results. As in
my
code, when N is odd, I have to use trial and error to figure out
which
one to use.

For example,

for N =6:

fftshift([-3:2]) = [ 0 1 2 -3 -2 -1 ]
ifftshift([-3:2]) = [ 0 1 2 -3 -2 -1 ]

whereas for N = 5

fftshift([-2:2]) = [ 1 2 -2 -1 0 ]
ifftshift([-2:2]) = [ 0 1 2 -2 -1 ]


It may appear trivial, but a one bin translation can cause an
imaginary
part with rapid ocscillations to appear when N is large.

>do I also need to do this with the variable?

No, I have just given you the variable indexing.

>I can evaluate the integral
>using residue calculus if I ignore the |k| terms and that should provide
> me with another check on my working.

???

Reread my post.

I have already found the poles without ignoring the abs(k).

> I had a look at you code and there were some commands which didn't run
> and that matlab.  Are you using an old version?

Yes. It is old as dirt.

>> ver
-------------------------------------------------------------------------------------
MATLAB Version 6.5.1.199709 (R13) Service Pack 1
MATLAB License Number: XXXXXX
Operating System: Microsoft Windows XP Version 5.1 (Build 2600:
Service Pack 3)
Java VM Version: Java 1.3.1_01 with Sun Microsystems Inc. Java
HotSpot(TM) Client VM
-------------------------------------------------------------------------------------

sortfigs is my own 3 line figure reordering function.

Hope this helps.

Greg

Subject: Inverse Fourier transforms

From: Mat Hunt

Date: 23 Jun, 2010 10:58:04

Message: 21 of 32

Hi Greg,

Still not really getting you. I have a very formal training in terms of Fourier analysis, it was a pure maths course I did in my third year at university.

I think of a Fourier transform as a map essentially which is defined by:

F(k)=int_{-infinity}^{+infinity}f(x)e^{-ikx}dx

The inverse transform being defined as:

f(x)=1/(2*pi)*\int_{-infinity}^{+infinity}F(k)e^{ikx}dk

In my work I have calculated what F(k), I now need to find what f(x) is.

You said that I need to do fftshift(ifft(F(k)) to get f, now I tried to read your post to figure out what I need to do but I couldn't quite get it. Can you explain again as if I were a simpleton :o)

The thing with the residue calculus was just a way of testing if I got the correct answer or not.

Mat

Subject: Inverse Fourier transforms

From: Greg Heath

Date: 23 Jun, 2010 17:38:27

Message: 22 of 32

On Jun 23, 6:58 am, "Mat Hunt" <hunt_...@hotmail.com> wrote:
> Hi Greg,
>
> Still not really getting you. I have a very formal training in terms of Fourier analysis, it was a pure maths course I did in my third year at university.
>
> I think of a Fourier transform as a map essentially which is defined by:
>
> F(k)=int_{-infinity}^{+infinity}f(x)e^{-ikx}dx
>
> The inverse transform being defined as:
>
> f(x)=1/(2*pi)*\int_{-infinity}^{+infinity}F(k)e^{ikx}dk
>
> In my work I have calculated what F(k), I now need to find what f(x) is.
>
> You said that I need to do fftshift(ifft(F(k)) to get f, now I tried to read your post to figure out what I need to do but I couldn't quite get it. Can you explain again as if I were a simpleton :o)

You said it. I didn't. ;>)

I repeat:

The MATLAB fft and ifft functions assume the nonnegative intervals x =
dx*(0:N-1) and k = dk*(0:N-1) corresponding
to the Fourier pair f(x) and F(k).

If f and F are considered periodic, fftshift yields the translation to
the "zero centered" bipolar intervals with
indexing

xb = dx * [-ceil((N-1)/2) : floor((N-1)/2)]
kb = dk * [-ceil((N-1)/2) : floor((N-1)/2)]

corresponding to the Fourier pair fb(xb) and Fb(kb).

ifftshift is the inverse of fftshift.

When N is even, fftshift and ifftshift yield
the same result.

However, when N is odd,

x = ifftshift(xb)
X = fft(x)
Xb = fftshift(X)

and to return

X = ifftshift(Xb)
x = ifft(X)
xb = fftshift(x)

The trick is to remember that fftshift will
"center" the waveform about the coordinate zero.

If in doubt type fftshift([-2:2]) into the command
line.

Hope this helps.

Greg

Subject: Inverse Fourier transforms

From: Greg Heath

Date: 23 Jun, 2010 17:53:19

Message: 23 of 32

On Jun 23, 1:38 pm, Greg Heath <he...@alumni.brown.edu> wrote:
> On Jun 23, 6:58 am, "Mat Hunt" <hunt_...@hotmail.com> wrote:
>
> > Hi Greg,
>
> > Still not really getting you. I have a very formal training in terms of Fourier analysis, it was a pure maths course I did in my third year at university.
>
> > I think of a Fourier transform as a map essentially which is defined by:
>
> > F(k)=int_{-infinity}^{+infinity}f(x)e^{-ikx}dx
>
> > The inverse transform being defined as:
>
> > f(x)=1/(2*pi)*\int_{-infinity}^{+infinity}F(k)e^{ikx}dk
>
> > In my work I have calculated what F(k), I now need to find what f(x) is.
>
> > You said  that I need to do fftshift(ifft(F(k)) to get f, now I tried to read your post to figure out what I need to do but I couldn't quite get it. Can you explain again as if I were a simpleton :o)
>
> You said it. I didn't. ;>)
>
> I repeat:
>
> The MATLAB fft and ifft functions assume the nonnegative intervals x =
> dx*(0:N-1) and k = dk*(0:N-1) corresponding
> to the Fourier pair f(x) and F(k).
>
> If f and F are considered periodic, fftshift yields the translation to
> the "zero centered" bipolar intervals with
> indexing
>
> xb = dx * [-ceil((N-1)/2) : floor((N-1)/2)]
> kb = dk * [-ceil((N-1)/2) : floor((N-1)/2)]
>
> corresponding to the Fourier pair fb(xb) and Fb(kb).
>
> ifftshift is the inverse of fftshift.
>
> When N is even, fftshift and ifftshift yield
> the same result.
>
> However, when N is odd,
>
> x  = ifftshift(xb)
> X  = fft(x)
> Xb = fftshift(X)
>
> and to return
>
> X  = ifftshift(Xb)
> x  = ifft(X)
> xb = fftshift(x)

WHOOPS! Screwed up the notation.
Correction;

 f = ifftshift(fb)
 F = fft(f)
 Fb = fftshift(F)

 and to return

 F = ifftshift(Fb)
 f = ifft(F)
 fb = fftshift(f)


> The trick is to remember that fftshift will
> "center" the waveform about the coordinate zero.
>
> If in doubt type fftshift([-2:2]) into the command
> line.
>
Hope this helps.

Greg

Subject: Inverse Fourier transforms

From: Mat Hunt

Date: 24 Jun, 2010 11:25:25

Message: 24 of 32

Hi Greg,

I had a play with it and to be honest, I can't get the thing to work. I don't understand the need for something like the fftshift function. I expected (naively) that I would give the point at which I would want the function evaluated at. Looking at the definition of the Fourier transform, they have the ability to do this, yet they choose not to for some reason.

I looked at what fftshift did to -2:2 and I have to say that I really didn't understand why it was doing that.

I looked at the maximum and minimum of the vector I got after doing ifft and it didn't match up to how I got doing it another method. I think I am going to keep with my old method of computing the inverse Fourier transform as it gives sensible. I wanted to compare my method with matlab's inbuilt method but after all this messing around, it still doesn't work.

I shall mention this to the developers as something that they might like to comsider.

Thanks for your time.

Mat

Subject: Inverse Fourier transforms

From: Mat Hunt

Date: 24 Jun, 2010 22:13:04

Message: 25 of 32

So can anyone give me a clear and simple methodology to go from a Fourier transformed function F(k), back to the original function at a point, x say? So this would be a single number? I would have thought the the people at matlab would have thought this one through.

Mat

Subject: Inverse Fourier transforms

From: dbd

Date: 24 Jun, 2010 23:21:02

Message: 26 of 32

On Jun 24, 3:13 pm, "Mat Hunt" <hunt_...@hotmail.com> wrote:
> So can anyone give me a clear and simple methodology to go from a Fourier transformed function F(k), back to the original function at a point, x say? So this would be a single number?  I would have thought the the people at matlab would have thought this one through.
>
> Mat

You persist in having continuous function expectations for discrete
finite implementations (fft, ifft). I think the people at Matlab
understand this There are symbolic and discrete implementations. With
data samples, you have picked the discrete. Abandon your symbolic
expectations. There are those who try to cure your problem with
fftshift and ifftshift, but frankly, until you understand that what
the fft and ifft functions give you are different from the infinite/
continuous theory of Fourier and how they are different, I doubt that
a review of fftshift and ifftshift syntax can have any useful meaning
to you.

Good Luck!

Dale B. Dalrymple

Subject: Inverse Fourier transforms

From: Mat Hunt

Date: 27 Jun, 2010 19:57:04

Message: 27 of 32

Okay, all I am asking is if I give a function F=F(k) which is the Fourier transform of a function f=f(x) and another point x, say then I have a number as an output which is f(x). Can the ifft function do this?

Mat

Subject: Inverse Fourier transforms

From: Greg Heath

Date: 28 Jun, 2010 06:24:59

Message: 28 of 32

On Jun 27, 3:57 pm, "Mat Hunt" <hunt_...@hotmail.com> wrote:
> Okay, all I am asking is if I give a function F=F(k) which is the Fourier transform of a function f=f(x) and another point x, say then I have a number as an output which is f(x). Can the ifft function do this?
>
> Mat

Yes.

However, (once again):

The fft and ifft functions assume that the waveforms are defined over

x = dx*(0:N-1); % dx = 1/Fs , Fs is the spatial sampling frequency

and

k = dk*(0:N-1); df = Fs/N

The defining summation equations in the documentation of fft and ifft
are periodic with period N.

Therefore if you start out with k data that are defined over
different intervals, e.g.,

kb = dk* [ -ceil((N-1)/2) : floor((N-1)/2 ]

you have to assume the function is periodic outside
of the given interval and find that portion that lies between
0 and dk*(N-1).

Hope this helps.

Greg

Subject: Inverse Fourier transforms

From: Mat Hunt

Date: 28 Jun, 2010 10:15:08

Message: 29 of 32

The problem is that my function isn't periodic. I am looking for a travelling wave solution which vanishes as x\rightarrow\pm\infty

Mat

Subject: Inverse Fourier transforms

From: Greg Heath

Date: 30 Jun, 2010 14:24:50

Message: 30 of 32

On Jun 28, 6:15 am, "Mat Hunt" <hunt_...@hotmail.com> wrote:
> The problem is that my function isn't periodic. I am looking for a travelling wave solution which vanishes as x\rightarrow\pm\infty
>
> Mat

99.9% of the waveforms that are used in DFT/FFT are
not periodic. However, when you use the DFT/FFT you
assume that, within a finite interval, the waveform
can be represented by a weighted sum of sines and
cosines of a fundamental frequency and it's harmonics.
Because of that, when that sum is evaluated outside
of the defining interval, it is periodic.

Therefore if you are using the DFT/FFT to represent
nonperiodic phenomena, you have to make that defining
interval sufficiently long.

For your physical problem it seems like there is
a system of linear partial differential equations
in x and t which you are trying to solve using dual
Fourier-Laplace transforms in space-time. I'd have
to know the complete mathematical specification of
the problem including initial and boundary conditions
before I can make any more relevant comments.

Hope this helps.

Greg

Subject: Inverse Fourier transforms

From: M Antola

Date: 13 Feb, 2013 15:49:12

Message: 31 of 32

Sorry for bringing up an old post, but this was a high rated link when googleing the problem and I didn't find a satisfactory solution. So I want to share the solution. Say I start with

>> syms k
>> g=@(k) k.^2.*exp(-k.^2)

g =

    @(k)k.^2.*exp(-k.^2)

This is the fourier space function, just a nontrivial example. Next, the analytic inverse:

>> ifourier(g(k))
 
ans =
 
((pi^(1/2)*exp(-x^2/4))/2 - (pi^(1/2)*x^2*exp(-x^2/4))/4)/(2*pi)
 
So we define this as f(x):

>> f=@(x) ((pi^(1/2).*exp(-x.^2/4))/2 - (pi^(1/2).*x.^2.*exp(-x.^2/4))/4)/(2*pi)

f =

    @(x)((pi^(1/2).*exp(-x.^2/4))/2-(pi^(1/2).*x.^2.*exp(-x.^2/4))/4)/(2*pi)

The main difficulty for me was to think of the optimal way to present the accuracy of the ifft sampling. You give ifft just a list of the value of g(k) at different values k_i. The vector k maps to a vector x automatically when we do ifft. But in a theoretical problem like this we probably want dk=dx, i.e. the sampling in dimensionless units to be the same in k and x space. So I did a little algebra to get the same sampling and range in both k and x space, resulting in the following piece of code:

deltatarget=0.1;
N0=1+2*pi/deltatarget^2;
N=2^nextpow2(N0);
L=sqrt(2*pi*(N-1));
x=linspace(-L/2,L/2,N);
k=x;

Now simply

fa=f(x); % the analytic inverse
fn=N/L*ifftshift(ifft(fftshift(g(k)),'symmetric')); % the numeric inverse

and using

plot(x,fa,'b',x,fn,'r')
xlabel('x')
ylabel('f(x), fn(x)')
axis([-5 5 -0.3 0.3])

we see the analytic and numeric inverse transforms are on top of each other.

Let me know if there is some way to improve.

Subject: Inverse Fourier transforms

From: M Antola

Date: 15 Feb, 2013 13:06:06

Message: 32 of 32

"M Antola" <matti.antola@helsinki.fi> wrote in message <kfgclo$4ot$1@newscl01ah.mathworks.com>...
> Sorry for bringing up an old post, but this was a high rated link when googleing the problem and I didn't find a satisfactory solution. So I want to share the solution.

I had a small mistake (worth of dx or dk) in my solution, so for completeness heres the fixed script. Also the algorithm is most accurate when the numeric x's start from minus Lx/2 (Lx=period) and end at Lx/2-dx, not the other way around, or the symmetric case.

f=@(x) exp(-x.^2)
g=@(w) pi^(1/2)*exp(-w.^2/4)

dk=0.01;
N0=1000;
Nx=2^nextpow2(N0);
dx=2*pi/(Nx*dk)
Lx=Nx*dx;
Lk=Nx*dk;

kn=linspace(-Lk/2,Lk/2-dk,Nx);
xn=linspace(-Lx/2,Lx/2-dx,Nx);

fa=f(xn);
fn=1/dx*ifftshift(ifft(fftshift(g(kn)),'symmetric'));

plot(xn,fa,'b',xn,fn,'r')
xlabel('x')
ylabel('f(x)')
axis([-5 5 -0.1 1.1])

sum((fn-fa).^2)/sum(fa.^2) % order e^-7 for these values.

Tags for this Thread

What are tags?

A tag is like a keyword or category label associated with each thread. Tags make it easier for you to find threads of interest.

Anyone can tag a thread. Tags are public and visible to everyone.

Contact us