Inverse fft doesn't match analytical formulation

I need to convert a function y=1/w^2 (y=0 at w=0) from frequency domain to real domain using inverse fourier transform.

The analytical solution of it is y=|x|/2. However, when I tried to use ifft, I cannot get the correct results.

I then manually coded a discrete inverse fourier transform and was able to reproduce the analytical formulation.

My own DIFT is slow so I still need the ifft. could anyone help identify the issue in my usage of ifft? Many thanks for your help!

Below are my code: (I didn't scale the results)

ffhandle = @(w) arrayfun(@(wi) (1 / wi^2) * (wi ~= 0) + (wi == 0) * 0, w); % define the function
%%% Total number of points in frequency domain (must be even)
N= 2^12;
%%% Computing length of reconstruction
L = 10;
%%% Delta Width for recovered data
W = L/(2*N);

m = ceil(L/W);
%%% Range of Frequency for DFT method
w0 = 31;
w1 = (-w0/2:w0/N:w0/2-w0/N);
fhat1 = feval(ffhandle,w1);
fhat1(N/2+1)=0; %remove singularity at zero

x = (-m*W-W/2:W:-W/2);
f = zeros(size(x));
dw = abs(w1(2)-w1(1));

for n = 1:length(w1)
    f = f + (dw/(2*pi))*fhat1(n)*exp(1i*w1(n)*x);     % this is my own ifft.
end
f2 = ifft(ifftshift(fhat1));   % this is matlab's ifft.
figure
plot(real(f))
hold on
plot(real(f2))
hold off

[Paul]

Hi shihao wang,

The code as provided has an error. Please provide the corrected code so we can replicate your results.

ffhandle = @(w) arrayfun(@(wi) (1 / wi^2) * (wi ~= 0) + (wi == 0) * 0, w); % define the function
%%% Total number of points in frequency domain (must be even)
N= 2^12;
%%% Computing length of reconstruction
L = 10;
%%% Delta Width for recovered data
W = L/(2*N);
m = ceil(L/W);
%%% Range of Frequency for DFT method
w0 = 31;
w1 = (-w0/2:w0/N:w0/2-w0/N);
fhat1 = feval(ffhandle,w1);
fhat1(N/2+1)=0; %remove singularity at zero
x = (-m*W-W/2:W:-W/2);
f = zeros(size(x));
dw = abs(w1(2)-w1(1));
for n = 1:length(w1)
    f = f + (dw/(2*pi))*fhat1(n)*exp(i*w1(n)*x);     % this is my own ifft.
end

Error using *
Incorrect dimensions for matrix multiplication. Check that the number of columns in the first matrix matches the number of rows in the second matrix. To operate on each element of the matrix individually, use TIMES (.*) for elementwise multiplication.

f2 = ifft(ifftshift(fhat1));   % this is matlab's ifft.
figure
plot(real(f))
hold on
plot(real(f2))
hold off

[shihao wang]

Hi Paul,

Thank for your prompt response and help on checking the code.

I have tried at my end using matlab 2022b and matlab online but didn't get the error message.

Actually for this line, f = f + (dw/(2*pi))*fhat1(n)*exp(i*w1(n)*x);

this is the discrete fourier transform and fhat1 stores the function values in the freqeuncy domain. w1 is the frequency array and x is the array of coordinates in the real domain. here only x is used for array-based calculation.

I suspect it maybe related to the version that f needs to be tranposed?

Could you try replace the line with f' = f' + (dw/(2*pi))*fhat1(n)*exp(i*w1(n)*x);? Please let me know if the error message persists. many thanks!

[Paul]

Well, this is very strange. For some reason i is being seen as an empty variable, instead of the the function that returns sqrt(-1).

which i

i is a variable.

whos i

  Name      Size            Bytes  Class     Attributes

  i         0x0                 0  double              

i

i =

     []

We can get around this for now by using 1i in the expression for f

ffhandle = @(w) arrayfun(@(wi) (1 / wi^2) * (wi ~= 0) + (wi == 0) * 0, w); % define the function
%%% Total number of points in frequency domain (must be even)
N= 2^12;
%%% Computing length of reconstruction
L = 10;
%%% Delta Width for recovered data
W = L/(2*N);
m = ceil(L/W);
%%% Range of Frequency for DFT method
w0 = 31;
w1 = (-w0/2:w0/N:w0/2-w0/N);
fhat1 = feval(ffhandle,w1);
fhat1(N/2+1)=0; %remove singularity at zero
x = (-m*W-W/2:W:-W/2);
f = zeros(size(x));
dw = abs(w1(2)-w1(1));
for n = 1:length(w1)
    f = f + (dw/(2*pi))*fhat1(n)*exp(1i*w1(n)*x);     % this is my own ifft.
end
f2 = ifft(ifftshift(fhat1));   % this is matlab's ifft.
figure
plot(real(f))
hold on
plot(real(f2))
hold off

Is that the plot you see on 2022b and Matlab Online?

[shihao wang]

Hi Paul,

Yes this is exactly what I get.

The blue line (from f) is closer to the correct answer but the built in ifft (red line) cannot get the same results.

It would be great if you can check the difference. Much appreciated!

Answer 1

[Torsten]

syms w
F = 1/w^2;
ifourier(F); disp(char(ans))

-(x*sign(x))/2

[shihao wang]

Hi Torsten,

Thanks for your help! The sysbolic operation is definetely a great way to crack the problem.

However, in my other problems I do need numerical ifft as the functions have no close-form expression.

Could you also help look at the implementation of my ifft? I think the difference lies in the treatment of frequency.

Thanks again and have a great day!

Answer 2

[Paul]

Let's take things one step at a time and see where they lead. If we have Continous-Time Fourier Transform (CTFT)

F(w) = 1/w^2

then, as @Torsten showed,

f(x) = -abs(x)/2.

So that should basically be our target function for reconstruction.

However, we have that pesky problem of what to do at w = 0, which you've decided to specify F(0) = 0.

In your code, there really is no need to define fhandle that way with the subsequent call to feval. Just define a vectorized, anonymous function

F = @(w) (1 ./ w.^2) .* (w ~= 0) + (w == 0) .* 0;

But that still doesn't work for F(0)

F(0)

ans = NaN

The reason that returned NaN is because

(1/0)*(0 ~= 0)

ans = NaN

What we need to do is turn that NaN into a zero (per your specification), which we can do using fillmissing, and we can shorten the expression because the (w == 0)*0 doesn't contribute

F = @(w) fillmissing( (1 ./ w.^2) .* (w ~= 0),'constant',0);

Let's check

w = [-0.7,0,0.3];
F(w)

ans = 1×3

2.0408 0 11.1111

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

1./w.^2

ans = 1×3

2.0408 Inf 11.1111

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

Continuing

w0 = 31;
N = 2^12;
%w1 = (-w0/2:w0/N:w0/2-w0/N);w1([1,end])
n = ((-N/2):(N/2-1));
w1 = n/N*w0;
fhat1 = F(w1);

At this point, we are going to define the sample points for reconstruction of f(x)

L = 10;
W = L/(2*N);
m = ceil(L/W);
x = (-m*W-W/2:W:-W/2);

I'm not quite sure what the intent was in defining x this way, let's make sure it's correct in terms of number of samples and their values. Note that x has more points than w1.

numel(x)

ans = 8193

x([1 end])

ans = 1×2

-10.0006 -0.0006

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

Should x be assymetric like that?

Continuing, it looks like the next section is the a numerical approximation to the inverse CTFT integral.

% this line is vectorized trapezoidal integration
%f = 1/2/pi*trapz(w1,fhat1.*exp(1j*w1.*x(:)),2);
% rectangular integration with a loop
f = zeros(size(x));
dw = abs(w1(2)-w1(1));
for ii = 1:length(w1)
    f = f + (dw/(2*pi))*fhat1(ii)*exp(1i*w1(ii)*x);     % this is my own ifft.
end

Need to plot against x

figure
plot(subplot(211),x,real(f)),grid
plot(subplot(212),x,imag(f))

Stopping here for now ...

[shihao wang]

Hi Paul,

Thanks for the step by step clarification.

For x, I intend to test an arbitrary range, because I want to check the ifft on any real domain. As such, x is no necessarily symmetrical in my case.

The slop of the real f plot is 0.5, but the built-in ifft is not. This is what I want to figure out.

Plz let me know if you need any further clarifications.

Thanks!

[Walter Roberson]

F = @(w) (1 ./ w.^2) .* (w ~= 0) + (w == 0) .* 0;

can be fixed to

F = @(w) 1./(w.^2 + (w ==  0)) - (w == 0)

F = function_handle with value:
    @(w)1./(w.^2+(w==0))-(w==0)

X = (-5:5)/10

X = 1×11

-0.5000 -0.4000 -0.3000 -0.2000 -0.1000 0 0.1000 0.2000 0.3000 0.4000 0.5000

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

F(X)

ans = 1×11

4.0000 6.2500 11.1111 25.0000 100.0000 0 100.0000 25.0000 11.1111 6.2500 4.0000

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

[Paul]

Continuing ...

F = @(w) fillmissing( (1 ./ w.^2) .* (w ~= 0),'constant',0);
w0 = 31;
N = 2^12;
%w1 = (-w0/2:w0/N:w0/2-w0/N);w1([1,end])
n = ((-N/2):(N/2-1));
w1 = n/N*w0;
fhat1 = F(w1);
L = 10;
W = L/(2*N);
m = ceil(L/W);

Define a symmetric x so we have an easier visual comparison to the target function -abs(x)/2

%x = (-m*W-W/2:W:-W/2);
x = linspace(-10,10,500);

% this line is vectorized trapezoidal integration
%f = 1/2/pi*trapz(w1,fhat1.*exp(1j*w1.*x(:)),2);
% rectangular integration with a loop
f = zeros(size(x));
dw = w1(2)-w1(1);
for ii = 1:length(w1)
    f = f + (dw/(2*pi))*fhat1(ii)*exp(1i*w1(ii)*x);     % this is my own ifft.
end

figure
plot(subplot(211),x,real(f)),grid
plot(subplot(212),x,imag(f))

The real part looks pretty reasonable wrt to the shape of the target function, with the vertical offset likely due to how we are handling F(0).

To process fhat1 through the ifft, we have to ifftshift, then take the ifft, then fftshift to center the result around x = 0.

f2 = fftshift(ifft(ifftshift(fhat1),'symmetric'));

In order to plot against the independent variable we need to derive the sampling period

Xs = 2*pi/N/dw;
x = n*Xs;

Now plot, keeping in mind we have to scale the ifft output by the sampling period to match up to the iCTFT.

figure
plot(x,real(f2)/Xs),grid

At first glance, that looks different than the previous result. But let's get everything on the same scale

x = n*Xs;
f = zeros(size(x));
dw = w1(2)-w1(1);
for ii = 1:length(w1)
    f = f + (dw/(2*pi))*fhat1(ii)*exp(1i*w1(ii)*x);     % this is my own ifft.
end

figure
plot(x,real(f),x,f2/Xs,'.','MarkerIndices',1:40:numel(x)),grid
hold on
plot(x,-abs(x)/2+max(real(f))),legend('iCTFT','iDFT','Target')

The recontructed signals are basically identical. The reconstruction only well-approximates the target function for an interval around x = 0.

Of course, we'll never get a perfect reconstruction over the entire real line because of the windowing of F(w) (-31/2 < w < 31/2), the sampling of F(w), the situation at F(0), and the imperfect integration of the frequency domain samples.

Not sure if that's the result you're looking for. We can take a different approach if you want a set of frequency domain samples for which the ifft matches the target function over a symmetric interval abs(x) <= x0. But I don't really know what the ultimate goal is of this exercise.

Also, keep in mind that with this approach the reconstructed signal is actually periodic. Here, we stretch x to cover three periods.

x = n*Xs;
x = x*3;
f = zeros(size(x));
dw = w1(2)-w1(1);
for ii = 1:length(w1)
    f = f + (dw/(2*pi))*fhat1(ii)*exp(1i*w1(ii)*x);     % this is my own ifft.
end
figure
plot(x,real(f)),grid

[shihao wang]

Hi Paul,

Thanks for your correction and explanation!

Could you explain a bit about this line: Xs = 2*pi/N/dw; x = n*Xs; what is the purpuse of the scaling? Also what is the purpose of the double fftshift?

Also, I notice the results are periodic. is there a way to avoid that?

To further clarify, some background of this practice is added here:

I want to solve odes using Fourier transform. As a first trial, I started with

In fourier domain, the above ode is solved as

because the delta function is 1 in the Fourier domain.

Then apply the inverse transform we get

I want to use discrete fourier transform to solve this problem. This is because in my other problems there is no close-form expression in the real domain.

Many thanks for your help!

[Paul]

Suppose we have a continuous time signal x(t) with Fourier tranform X(w). The Fourier transform of the uniformly-impulse-sampled version of x(t) with sampling period Xs is X_p(w), which is the periodic extenstion of X(w) with period 2*pi/Xs. Now if we take N uniformly spaced samples of X_p (w) over one period, those samples will be spaced by dw = 2*pi/Xs/N.

One interpreration of fft/ifft is that the inputs to those functions are one period of a periodic signal and they each return one period of a periodic signal. By convention, the input to fft/ifft is the period of the signal that starts at 0. Here, we have fhat1 defined over -31/2 < w < 31/2 as one period of a periodic signal. So on input to ifft we have to ifftshift to get the single period from 0 < w < 31, which you had in your original code. Now the output of ifft is one period of the signal defined over 0 < x < N*xs, which is basically what you got, i.e., the period of the output over 0 < x < 900 (approximately), refer to the "three period" plot in my previous comment. So we have to fftshift the output of ifft to get to the "central" period, which is what we want.

I don't think the periodicity can be avoided once we go down the sampling path. Sampling in one domain, (frequency in this case), results in periodicity in the other domain (space in this case).

There is a sign error in your trial problem. It is true that Y(w) = -1/w^2, but the inverse transform of that is: abs(x)/2 (w/o a negative sign).

So you have closed-form expressions for the Fourier transform Y(w), but you can't find the closed form expression for y(t)? Can you give a simple example of such a Y(w)?

[Paul]

An approach that seems to work o.k. for this problem, within a constant, is to shift in the frequency domain to avoid sampling the frequency domain signal at the singularity.

F = @(w) -1./w.^2;  % based on the example differential equation above

n = 2^12;
N = 2*n+1;
w0 = 31;
w = (-n:n)*w0/N;
dw = w(2) - w(1);
dx = 2*pi/N/dw;
x = (-n:n)*dx;

Fshift = @(w) F(w - dw/2); % shift by half of a dw

Shift in frequency domain is multiplication by exp in the time domain

fhatshift = Fshift(w);
f = fftshift(ifft(ifftshift(fhatshift))).*exp(-1j*dw/2*x)/dx;

Imaginary part of f is small, throw it away

max(abs(imag(f)))

ans = 3.8550e-11

f = real(f);

figure
plot(x,f,x,abs(x)/2),grid

If we posit that the true answer should be 0 at x = 0, then we can offset f by

sum(fhatshift)/N/dx

ans = -415.1257

However, if we compare to the target fuction and zoom in very tight we do see that the reconstruction isn't perfect and has oscilations.

figure
plot(x,f-min(f)-abs(x)/2)
ylim([-0.0205362,-0.0205360])

Answer 3

[David Goodmanson]

Hi SW,

There is a sign error in your transform pair, which should be

1/w^2 <--> -|x|/2   not  1/w^2 <--> |x|/2

As has been noted, that transform pair is not well suited to the fft & ifft. That's because

o The fft represents a periodic function that repeats with a periodicity of the width of the window (as in Paul's example above). That periodicity is inevitable. If the function is, for example a cosine wave that fits exactly into the window, then the periodicity is natural. At other times, say with a parabola, it is forced.

o Any function can be made to repeat, but if f(x) at the left boundary of the window and f(x) at the right boundary of the window do not agree in both value and derivatives, that represents a discontinuity in f(x). The discontinuity can affect the fourier transform and obscure the result you are looking for, sometimes by a lot. In the present case the |x| values do agree on the right and left boundaries, but their first derivatives do not.

o If the both function and its transform fall off to negligible values at the edges of the window, such as the pair of gaussians e^(-x^2) and e^(-w^2), then fft & ifft will (usually) work to good precision. In the present case while 1/w^2 can be made to fall off reasonably well, obviously |x| does not.

Even if the fft & ifft didn't have those deficiencies, I don't think the technique of setting the value of 1/w^2 to 0 at w = 0 is a good way to go about it. For a function g(w) that approximates 1/w^2, you want the transform

Int{-inf,inf} g(w) e^(iwx) dw = |-x|/2.  

When x = 0, we have

Int{-inf,inf} g(w) dw = 0.

and this can only happen if the integrand g(w) has regions of both signs. Setting 1/w^2 to 0 at w = 0 leaves an integrand that is positive everywhere else on the w axis. So you cannot possibly attain an integral of 0 when x = 0. The removal idea is on the right track but method below does something similar in a more predictable way.

In the following, continuous functions are used as a model for what goes on, with the normalization

Int f(x) e^(-iwx) dx = g(w)            (1/2pi) Int g(w) e^(iwx) dw = f(x)       (3)

and all integrals are from -inf to inf.

First of all, the transform

Int 1/w^2 e^(-iwx) dx = -|x|/2

can't really be true because for x = 0 you are looking for a result of 0, but the integral of 1/w^2 is infinite. There has to be a limiting process involving functions f(x,eps) and g(k,eps) such that

     as eps-> 0, f(x,eps) ->  -|x|/2,  for fixed x                          (e1)
and  as eps-> 0, g(w,eps) ->   1/w^2,  for fixed w                          (e2)
        
examples  (e1) (-1/2) sqrt(x^2 + eps^2)
          (e2) 1/(w^2 + eps^2)

Define h(x) as the heaviside function and let

          f(x) = h(x) e^(-eps x)

which is a falling expontential for x>0. The transform pair is

h(x) e^(-eps x) <-->  1/(eps + iw)

and for the rising exponential for x<0, the transform pair is

h(-x) e^(eps x) <-->  1/(eps - iw) 

(No actual work to get the second one since for a given transform you can get another transform with x -> -x & w -> -w).

Adding these we arrive at the pair

e^(-eps |x|) <-->  2eps/(w^2 + eps^2) 

If both sides are divided by 2eps, the right hand side goes as (e2) above, but the left hand side doesn't work. Something similar to this has been used on this thread, but the results show functions that do not equal 0 when x = 0.

The situation can be fixed by taking the derivative of this transform pair w/r/t eps, The result is the pair

          f(x)                             g(w)

    (-|x|/2) e^(-eps |x|) <-->  (w^2 - eps^2)/(w^2 + eps^2)^2 

which satisfies both (e1) and (e2). And it can be checked with ifft, remembering the limitations of fft & ifft.

The first plot shows an expanded view of g(w), which does integrate to zero. The positive parts fall off like w^-2, and the large negative spike is reminiscent of what you were doing when setting the function equal to 0 at w = 0. The domain of w is -1e4 to 1e4, so this plot is a closeup indeed.

The second plot is f(x), which does behave as -|x|/2 at the origin as you can see. The function tries, but before too long it falls away from the -|x|/2 line.

% w domain function
delw = .1;                 % w grid spacing
n = 1e5;
w = (-n:n)*delw;
N = 2*n+1;                 % number of fft points 
delx = 2*pi/(delw*N);      % x grid spacing;   fft golden rule is delx*delw = 2*pi/N
eps = 1/3;
g = (w.^2-eps^2)./(w.^2+eps^2).^2;
fig(1)
plot(w,g,'o-')
title('g(w)  closeup')
xlim([-10 10])

% x domain function
x = (-n:n)*delx;
y = fftshift(ifft(ifftshift(g)))*N*(1/(2*pi))*delw;
% factor of N since ifft automatically multiplies by 1/N and we are doing
%        things differently.
% factor of 1/2pi from normalization (3) in the text.
% factor of delw to approxmate the ifft sum by an integral.
% ifftshift puts the point corresponding to w=0 at the bottom of the array
% where ifft needs it to be; fftshift puts the x = 0 point into the middle
% of the array for display purposes
fig(2)
subplot(2,1,1)
plot(x,y) 
subplot(2,1,2)
plot(x,y,x,-abs(x/2))
legend('y(x)', '- abs(x)/2')
title('y(x)   closeup')
xlim([-.03 .03])