Inverse fft doesn't match analytical formulation
Show older comments
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
4 Comments
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
f2 = ifft(ifftshift(fhat1)); % this is matlab's ifft.
figure
plot(real(f))
hold on
plot(real(f2))
hold off
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
whos 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
on 2 Dec 2024
Answers (3)
Torsten
on 2 Dec 2024
Edited: Walter Roberson
on 3 Dec 2024
syms w
F = 1/w^2;
ifourier(F); disp(char(ans))
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
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)
The reason that returned NaN is because
(1/0)*(0 ~= 0)
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)
1./w.^2
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)
x([1 end])
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 ...
6 Comments
shihao
on 4 Dec 2024
F = @(w) (1 ./ w.^2) .* (w ~= 0) + (w == 0) .* 0;
can be fixed to
F = @(w) 1./(w.^2 + (w == 0)) - (w == 0)
X = (-5:5)/10
F(X)
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
on 5 Dec 2024
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
on 6 Dec 2024
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)?
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)))
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
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])
David Goodmanson
on 6 Dec 2024
Edited: David Goodmanson
on 6 Dec 2024
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])
Categories
Find more on Frequency Transformations in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!
