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

# Thread Subject: Inverse Fourier transforms

Subject: Inverse Fourier transforms

From: Mat Hunt

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

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

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

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

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

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

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

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

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

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

### 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
really trying to do.

Subject: Inverse Fourier transforms

From: Greg Heath

### 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 ### 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 ### 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 ### 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 ### 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 ### 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 ### 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 ### 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 ### 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 ### 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 ### 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 ### 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 ### 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 ### 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" 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 ### 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 ### 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 ### 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 ### 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 ### 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 ### 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.