Dft and continuous Fourier transform

[eoghan]

Hi there!
I need to calculate the Fourier transform of a continuous function in C++. To do this I need to use the Dft, but what is the relation between the Dft and the continuous Fourier transform? I mean, how can I get the continuous Fourier transform from the Dft?

 

[marcusl]

You have options.
1. If your function f(t) is continuous and closed-form, you might find its exact FT in a math book such as the CRC Standard Math Tables and Formulae. The SW packages Mathematica and Maple can also take continuous FTs.

2. If you want to take the DFT F(w_k) instead, make sure you bandlimit your function f(t) and sample it properly according to Shannon's Sampling Theorem. You can interpolate between the frequencies w_k using the appropriate sin(Nw)/sinw Dirichlet function.

3. If you just want the spectrum in a small range of frequencies, perhaps at closer spacing than you get with the DFT, use a chirp z transform with the radius factor set to 1.

You can find more info on the latter two approaches in any DSP text.

 

[eoghan]

Uhm.. my problem is this: I have to calculate the Fourier transform of a gaussian centered in x=5, i.e. I have to calculate:
$$\hat x(f)=\int_{-\infty}^{+\infty}e^{-(x-5)^2}e^{-ifx}dx$$
To do this I sample the function in a window large T=50 with N=100 sampling points, i.e., I sample the function from x=0 to x=50 with steps of T/N=0.5. In this way I get a discrete sequence of numbers, so that I can apply the DFT. After the Dft is applied, I get another sequence of numbers, let's call it $$x_{dft}$$. Now, is it true that the continuous Fourier transform is:
$$\hat x=x_{dft}(\frac{2\pi m}{T})$$
for m=0,..., N-1?

 

[Petr Mugver]

You can calculate the Fourier transform of the gaussian analitically... is there a particular reason for which you want the computer to do it for you?

 

[eoghan]

Yes, I have to solve the Schroedinger equation numerically using the split operator technique. In order to test the program I use a gaussian so I know exactly what the result should be. My goal is to study numerically the behavior of an electron in an harmonic oscillating e.m. field and in that case the Fourier transforms cannot be calculated analitically.

 

[marcusl]

The DFT is given by

$$F_k = \displaystyle\sum_{n=0}^{N-1}f_n \, exp{\left(-\frac{i2\pi nk}{N}\right)}$$

You have some issues to consider regarding sampling your function. First of all, a Gaussian is not bandlimited (it has energy extending to infinite frequencies). You have two choices: filter your time sequence so it's strictly bandlimited, or ignore it and accept that you'll have a bit of aliasing (error in your spectrum). Either way your result will be close but not exact.

Second, your choice of t=0 to 50 chops off the negative time portion of your waveform. This is a problem because the DFT builds into the spectrum the assumption that the signal is periodic. Your periodic signal will have a discontinuity at t=0, 25 that introduces high frequency ringing into your spectrum. Choose a time window that's symmetric about t=5 (the Gaussian peak) to avoid this. You can still index the samples as n=0 to 99.

 

[eoghan]

Ok, choosing a window that's symmetric about t=5 (I sample from t=-20 to t=25) I get the right Fourier transform. The problem now is when I transform back: I perform the inverse dft on the transformed sequence which is already discrete. I should obtain the initial gaussian centered at t=5, but I get a strange function (a kind of U centered at 0). This is the code:

	double x;
	int i;
	
	for(i=0, x=-20; i<N; i++, x+=(T/N))	//Sample the Gaussian with time-step T/N
	{	
		in2[i][0]	=	20*exp(-(x-5)*(x-5));          //Real part
		in2[i][1]	=	0;                            //Imaginary part
	}
	fftw_execute(fwd);   //Calculate the direct dft and store the result in out2

	for(i=0, x=0; i<N; i++, x+=T/N)	//Copy the result from out2 to in1
	{	
		in1[i][0]	=	out2[i][0];   //Real part
		in1[i][1]	=	out2[i][1];   //Imaginary part
	}
	fftw_execute(rew);              //Execute the inverse Dft
		
	
   //Draw the result with the function glVertex3f(x,y,z)
	glBegin(GL_LINE_STRIP);
		for(i=N/2; i<N; i++) 	glVertex3f((i-N)*2.*PI/T, 1./N*sqrt(out2[i][0]*out2[i][0]+out2[i][1]*out2[i][1]), 0);
		for(i=0; i<=N/2; i++) 	glVertex3f(i*2.*PI/T, 1./N*sqrt(out2[i][0]*out2[i][0]+out2[i][1]*out2[i][1]),0);
	glEnd();

Can you please check it out? Thank you very much!

 

[marcusl]

It's a well-known property of the FFT algorithm that the output spectrum is circularly rotated so as to extend from 0 to 2*pi/T_samp instead of from -pi/T_samp to +pi/T_samp. If you're a solid state guy, it's equivalent to defining the 1st Brillouin zone as 0 to 2*k_max instead of -k_max to +k_max. The information is all there and the same either way. The solution is to rotate your output data set halfway around. (Matlab has a command dedicated to this called fftshift.)

It sounds like your C++ routine shifts only on the inverse FFT. The solution (rotate the data) is the same.

 

[eoghan]

Uhm.. it's more complicated than I thought... do you know a good book (or a web site) where all this stuff is explained?

 

[marcusl]

It's more complicated to explain than to do. Look at your data set as circular, that is with the beginning and end connected together, and just rotate halfway. Your "U" will turn back into a Gaussian.

I learned years ago out of Oppenheim and Schafer, Digital Signal Processing. The latest incarnation of that text is called Discrete-Time Signal Processing. I don't have personal experience with other books but many people like Brigham, The Fast Fourier Transform. Most books on DSP will cover this material as well.

Here are the first two links that came up when I googled "fftshift":

http://www.ele.uri.edu/~hansenj/projects/ele436/fft.pdf"
http://www.dsprelated.com/showmessage/20790/1.php"

 

[eoghan]

Ok, I rotated the set of data and now I get a Gaussian, but it's not centered at x=5, but as far as I have to calculate this is not important. Anyway... thank you for your help! I'll have a look at the books!

 

[marcusl]

It's centered because you chose an original data set that was symmetric. The DFT has no knowledge of absolute time, only relative time starting from the first sample.