1
$\begingroup$

I am following my optics textbook (Optics, by Eugene Hecht), throughout which are given various exact analytical results for the diffraction patterns that result from light passing through differently shaped apertures. These diffraction patterns are often calculated by taking continuous Fourier transforms. Ultimately, I am trying to get to the stage where I can use a FFT to numerically calculate these patterns when an analytical solution is not possible. Therefore, as a first sanity-check step I would like to check that my FFT results match up with the analytical Fourier transforms in some simple cases including the correct scaling factors.

Unfortunately, I have failed at the first hurdle, in the simplest possible case I can think of!

The book gives the following definition for the continuous Fourier transform pair:

\begin{align} f(x) &= \frac{1}{2\pi} \int_{-\infty}^{+\infty} F(k)\; e^{-ikx}dk\hspace{1cm} &\text{Inverse Transform} \\ F(k) &= \int_{-\infty}^{+\infty} f(x)\; e^{ikx}dx\hspace{1cm} &\text{Forward Transform} \end{align}

The book then goes on to show some examples of functions along with their transforms - in particular, the two simple cases shown below. In Figure (a) it is shown that the spectrum of a constant $f(x) = A$ has a peak located at $k=0$. The peak is shown as an arrow, with a label of $2\pi A$:

$\hskip1in$ enter image description here

I tried to arrive at this by proceeding as follows:

\begin{align} F(k) &= \int_{-\infty}^{+\infty} f(x)\; e^{ikx}dx \\ &= \int_{-\infty}^{+\infty} A\; e^{ikx}dx \\ &= A\int_{-\infty}^{+\infty}\cos(kx)\;dx\; +\; iA\int_{-\infty}^{+\infty}\sin(kx)\;dx \end{align}

Now, because the cosine and sine waves extend from $\pm\infty$, I know that their areas will average out and the integrals will evaluate to zero EXCEPT for the case when $k=0$. In this special case, the sine term will still evaluate to zero (because $\sin(0)=0$), however the cosine term will evaluate to $\infty$, because the area under $\cos(0)=1$ will be infinite. This behaviour is captured by the Dirac delta function: $$\delta(k) = \left\{ \begin{array}{1 1} +\infty & \quad k = 0\\ 0 & \quad k\neq 0 \end{array} \right.$$ Therefore, I would say that the transform should evaluate to: $$F(k) = A\;\delta(k)$$ and when $k=0$ this should take a value of $+\infty$. I don't understand (a) why the value is finite, and (b) why there is a factor of $2\pi$ in the book's graph.

I also then tried to calculate numerically the value of the spectrum coefficient using Matlab's FFT function, which implements the DFT, using a value of $A=3$. According to the textbook, I guess I should see a value of $2\pi\cdot 3=18.8$.

num_points = 10;
dx = 0.5;                     % Sampling interval
x = 0:dx:(dx*num_points-dx);  % Generate position vector
A = 3;                        % Value of constant function
f = A*ones(num_points,1);     % Input function f(x) = A
F = fft(f)*dx;                % Scale FFT

disp(['DC Component = ',num2str(F(1))])

I understand there are some disagreements between people for right way to scale the FFT output (see here for example). Since the DFT is defined as $$Y_k = \sum^{N-1}_{n=0}y_n\; e^{-2\pi i k n/N}$$

it seems to make sense to me that this is a summing approximation to the continuous integral version of the FT above, and so should be scaled by $dx$ in order to give results that can be compared directly - so that is what I have done in the Matlab snippet. Unfortunately the script shows a value for the DC peak of 15 (and depends on the number of samples, which must be wrong).

Can anyone help me resolve these conflicts? As I say, I am looking to be able to calculate numerical FFTs and have them match the continuous analytical versions. I would like to move on to try a more complex function that can still be done analytically (such as a pure cosine, or a Gaussian), but if I can't do the simple constant value I don't think I should move on just yet!


Similar SE questions for reference:

Fourier Transform and Delta Function

About Fourier transform of periodic signal

Analytical Fourier transform vs FFT of functions in Matlab

how fft points should be rescaled to get the same results as the analytical solution?

Matlab FFT for gaussian function

$\endgroup$
1
$\begingroup$

Welcome! Regarding your questions:

a) Its value is not finite. $A$ (or $2\pi A$, which is the correct answer) is not its value, it's the delta function's "area". $\delta(k)$ is not a proper function, it's a distribution or a generalized function. When we sketch such functions, we denote their area as a scaling factor. Indeed the actual value of the delta function is $+\infty$ at $k=0$.

b) You recognize that the forward Transform is somewhat problematic. Try the inverse one to see that the FT pair holds: $$f(x) = \frac{1}{2\pi}\int_{-\infty}^{+\infty}F(k)e^{-jkx}dk = \frac{1}{2\pi}2\pi A\int_{-\infty}^{+\infty}\delta(k)e^{-jkx}dk= Ae^{-jkx}\Big|_{k=0} = A$$

So when you use MATLAB to approximate this pair, you can never get the $A=3$ number. I would use another piece of code to approximate the FT of the constant, which gives better control on scaling etc.

 dx = 0.01;              % x-step
 x = -10:dx:10;          % x-axis
 A = 3;                  % constant value
 f = A*ones(size(x));    % just an approximation of the constant
 dk = 0.01;              % k-step
 k = -10:dk:10;          % k-axis
 F = dx*f*exp(j*k'*x).'; % Riemann approximation of CTFT
 plot(k, F);             % plot

plot

As you can see from the plot, you can only approximate constant $A$ by some time-limited rectangular pulse. In the $k-$domain you get something that approaches the behavior of the delta function more and more as you increase the samples of your x-axis and the duration of your rectangular pulse.

You'll never get a perfect dirac delta because you would need infinitely many samples of your $k$ and $x$ axes. What you really get in the transform is the so-called sinc function that approaches the delta function under the notion of limit.

| improve this answer | |
$\endgroup$
  • $\begingroup$ Thanks, that is a little clearer. I still have a couple of questions though. (1) I totally see that performing the inverse transform works and gives back the original $f = A$ function. But I don't see why the $2\pi$ appears on the forward transform. Did I make a mistake in my working, or is the definition from the book wrong? I shouldn't need to perform the inverse transform, and then manually add afterwards a fudge factor of $2\pi$. $\endgroup$ – teeeeee Mar 6 at 10:54
  • $\begingroup$ (2) I don't understand why I should have to use the Riemann sum - is this something peculiar to the fact that I am trying to transform a limiting case. Would this problem go away if, for example, I was looking at a Gaussian function. Then the FFT should be able to reproduce the result of the definition of the CFT given in the textbook? Thanks for your patience. $\endgroup$ – teeeeee Mar 6 at 10:56
  • $\begingroup$ (1) The $2\pi$ is due to the scaling factor $1/2\pi$ of the inverse Transform. You can see that in the Fourier Transform f(x) is synthesized by exponentials of amplitude $1/2\pi$ times $F(n\Delta k)\Delta k$, not just $F(n\Delta k)\Delta k$. (2) You don't have to use it, it's just my most convienient way to handle CTFTs. You can do similar things with FFT, but different software packages have FFT implementations that sometimes differ in a scaling factor. Writing your own CTFT gives you full control on what's going on. Riemann approximation is one way to do it. $\endgroup$ – GKH Mar 6 at 21:25
  • $\begingroup$ Please forgive my misunderstanding, but if I want to calculate a transform of a function and don't know what its transform will be apriori (as in the case above), then I would never have known to add an extra factor of $2\pi$. I mean, if I only had the definition of the forward transform to work with, then my method of calculating the transform would never lead to a factor of $2\pi$, and I would arrive at the wrong result. Can you point to excatly which step in my proof goes wrong? Thanks! $\endgroup$ – teeeeee Mar 9 at 12:48
  • $\begingroup$ I agree that indeed as you say $f(x)$ is synthesised from exponentials with a $2\pi$ factor, but $F(x)$ is not synthesised from $f(x)$ with any $2\pi$.. $\endgroup$ – teeeeee Mar 9 at 12:50

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.