Fourier transformation is a classical and powerful tool in mathematics. However, there are also some details need to be given special attention during the numerical computation.

Principle of Fourier transformation

Fourier transformation is to convert a signal in the time/space domains $sig(t),\ sig(x)$ to the domains of frequency or wave number $sig(f),\ sig(k)$ . The Fourier transformation is defined as [1]:

sig(f) = \int sig(t)e^{i\omega t}dt

The related inverse Fourier transformation is:

sig(t) = \frac{1}{2\pi}\int sig(\omega)e^{-i\omega t}df

Where the signal frequency $f$ has the relation with the angular frequency/velocity as $\omega = 2\pi f$ .

Discrete Fourier Transformation

In numerical computation, we must use discrete data, thus we must convert the analytic Fourier transformation to discrete transformation. Which is to say, we need to convert the integration to summation as [2]:

sig_{fk} = \sum_{n=0}^{N-1} sig_{tn}e^{2\pi ikn/N}

The related inverse discrete Fourier transformation is:

sig_{tn} = \frac{1}{N}\sum_{k=0}^{N-1} sig_{fk}e^{-2\pi i kn/N}

Fast Fourier Transformation (FFT)

Fast Fourier Transformation uses a highly efficient approach to calculate the discrete Fourier transformation. In many programming language, there are well written FFT module to use, such as matlab, python and fortran. However, during the calculation of Fourier strength, careful attention should be given. A standard example is given by matlab help doc as [3] shows. Among this code example, several key points shall be noticed.

average fft amplitude by dividing the length of fft window L

1 2	Y = fft(sig) P2 = abs(Y/L)

This is because the fft approach by default do not normalized the strength. The raw result is a summation, the longer the nfft window, the stronger the strength. Thus, we must normalized the frequency amplitude by dividing the nfft window length [4].

correct real half frequency amplitude

1
2
3

P2 = abs(Y)
P1 = P2(1:L/2+1);
P1(2:end-1) = 2*P1(2:end-1);

Since the positive and negative frequency are joined together in the half positive strength, we should multiply the non-zero frequency amplitude by 2. However, the zero frequency component is unique, thus we don’t need to multiply it [3].

We can review this problem with python code as:

import numpy
import pylab

pylab.close('all')
pi = numpy.pi

Fs = 1000
dt = 1.0/Fs
t = numpy.arange(1, 3, dt)
f1 = 50
f2 = 300
sig = 3*numpy.sin(2*pi*f1*t) + 2*numpy.sin(2*pi*f2*t) + 1

nfft = len(t)
sig_fft = numpy.fft.fft(sig)
sig_fft_full = sig_fft/nfft
f = numpy.linspace(0, Fs/2, numpy.int(nfft/2))
f_just = numpy.append(f, numpy.linspace(-Fs/2, 0, numpy.int(nfft/2)))

# reconstruct the signal with FFT data
sig_t = numpy.fft.ifft(sig_fft)

# correct the half spectrum strength
sig_fft_Amp = numpy.abs(sig_fft)/nfft
sig_fft_real = sig_fft_Amp[0: numpy.int(nfft/2)+1]
f_real = numpy.linspace(0, Fs/2, numpy.int(nfft/2)+1)
sig_fft_real[1:numpy.int(nfft/2)] = 2*sig_fft_real[1:numpy.int(nfft/2)]

pylab.figure(figsize=(6, 6))
ax1=pylab.subplot(3,1,1)
pylab.plot(t, sig, label='sig(t) orignal')
pylab.xlabel('t (s)')
pylab.xlim([2, 2.1])
pylab.legend(loc='upper left', fontsize=10)
pylab.title('FFT and inverse iFFT')

ax2 = pylab.subplot(3,1,2)
pylab.plot(f_just, numpy.abs(sig_fft_full), '--', markersize=8, color='blue', label='sig(f) spectrum')
pylab.plot(f_real, sig_fft_real, 'o', color='red', label='real half spectrum')
pylab.legend(loc='upper left', fontsize=10)
pylab.xlabel('f (Hz)')
pylab.xlim([-Fs/2, Fs/2])

ax3 = pylab.subplot(3,1,3)
pylab.plot(t, sig, '--', color='blue', label='sig(t) original')
pylab.plot(t, sig_t, '-.', color='red', label='sig(t) reconstruct')
pylab.legend(loc='upper left', fontsize=10)
pylab.xlabel('t (s)')
pylab.xlim([2, 2.1])

pylab.tight_layout()
ax1.text(2.003, -2.2, '(a)', fontsize=18, color='black')
ax2.text(-Fs/2+30, 0.5, '(b)', fontsize=18, color='black')
ax3.text(2.003, -2.2, '(c)', fontsize=18, color='black')

Calculation result:

FFT, iFFT and amplitude correction

In this example, the test signal is: $sig(t)=3*sin(2*pi*50*t)+2*sin(2*pi*300*t)+1$ , this means that in the signal the strength of the 50 Hz component is 3, the amplitude of 300 Hz component is 2, the amplitude of 0 Hz component is 1. With the above code calculation, we can clearly see that after correction in the real half spectrum, the amplitude of the frequency component match to the given signal very well as figure (b) shows. Also the reconstructed signal matches well with the original signal. After the normalization, the unit of the spectrum density shall be the same as the original signal. eg, Bp(t) in Gs, Bp(f) also in Gs.

Reference

[1] 中科大高数教研室，高等数学导论，中科大出版社，合肥，2008，第三版。

[2] http://mathworld.wolfram.com/DiscreteFourierTransform.html

[3] https://www.mathworks.com/matlabcentral/answers/102727-why-does-the-example-in-fft-divide-the-fft-result-by-the-length-of-the-signal-and-multiply-the-magni

[4] https://in.mathworks.com/matlabcentral/answers/345736-why-fft-amplitude-changes-when-the-signal-s-length-changes