18.4.1.2 Algorithms (FFT)


A discrete Fourier transform (DFT) converts a signal in the time domain into its counterpart in frequency domain. Let (x_i) be a sequence of length N, then its DFT is the sequence (F_n) given by

F_n=\sum_{i=0}^{N-1}x_ie^{-\frac{2\pi j}{N}ni}

Origin uses the FFTW library to perform Fourier transform. With the transformed data, the amplitude, magnitude and power density can be computed by Origin.

FFTW

In FFTW, the computation of the transformed data is performed by an executor that is comprised of blocks of C code called "codelets". Each codelet specializes in one part of the transformation. With these codelets, the executor implements the Cooley-Turkey FFT algorithm, which factors the size of the input signal. By recursive factoring, the signal is broken into shorter parts. The results of the transforms of the short parts are multiplied; and finally the transform of the original signal is computed. More information on FFTW is available at http://fftw.org/.

Power density

By definition, power density or spectrum can be computed with the following equation:

P_{xx}(e^{j\omega })=\sum_{m=-\infty }^\infty r_{xx}(m)e^{-j\omega m}

where r_{xx}(m)\,\! is the auto-correlation function of the input signal.

However, we have a finite number of samples for the input signal and, therefore, computing the power spectrum with the definition is not possible as only certain methods can be used to estimate the power spectrum. The method used in Origin is the Periodogram, which estimates the power from the amplitude of the Fourier transformed data. While it is generally accepted that the squared amplitude is in proportion to the amplitude of power spectrum, various conventions exist for describing the normalization of power spectrum in each domain. The three descriptions used by Origin are "mean squared amplitude"(MSA), "sum squared amplitude"(SSA) and "time-integral squared amplitude"(TISA). They can be expressed as follows:

Power Density(two-sided)=\begin{cases}\frac{{Re}^2+{Im}^2}{n^2},for MSA\\\frac{{Re}^2+{Im}^2}n,for SSA\\\frac{\Delta t({Re}^2+{Im}^2)}n,for TISA\end{cases}

where Re\,\! and Im\,\! are the real and imaginary parts of the transform data; n\,\! is the length of the input sequence; \Delta t\,\! is the sampling interval.

Power spectrum can be either one-sided (single-sided) or two-sided (double-sided), depending on whether Two-Sided (2) or One-Sided (1) is selected for Spectrum Type (st). To compute the one-sided power density, it is necessary to compute the two-sided power density first. The result is then converted to the one-sided power using the following equations:

P_s(i)=P_d(i),i=0\,\!

P_s(i)=2P_d(i),i=1,2,\cdots \frac n2-1

where P_s(i)\,\! is the one-sided power spectrum and P_d(i)\,\! is the two-sided power spectrum.

If a window function is applied, the power result will be multiplied by a factor for compensation which is defined by :

N/{\sum_{n=0}^{N-1}w(n)^2}, where w(n) is the window function defined below.

More results

Origin can calculate the magnitude, phase and amplitude of the transformed data. Let Re\,\! and Im\,\! be the real and imaginary parts of the transform data, and let n\,\! be the size of the input signal. Use \Delta t\,\! to represent the sampling interval. Suppose the norma variable is set to 0 (normalization is not used). More outputs are calculated with the following formulas:

Spectrum Type
is Two-sided
(i=1-n/2 ~ n/2)
Spectrum Type
is One-sided
(i=0 ~ n/2)

Phase

\arctan (\frac{Im}{Re})\,

Magnitude

\sqrt{Re^2+Im^2}\,

Amplitude

\sqrt{Re^2+Im^2}/n\,

\sqrt{Re^2+Im^2}/n, i=0\mbox{ or }i=n/2\,
2*\sqrt{Re^2+Im^2}/n, \mbox{ otherwise }\,

dB

20log(Amplitude)\,

Normalized Amplitude in dB

dB-max(dB)\,

RMS amplitude

\frac{\sqrt{2}}2Amplitude\,

Normalization

The above computations are actually based on the assumption that the norma vairable is set to false. If this variable is set to true, the complex, real, imaginary, magnitude and square magnitude results will be normalized. Note that phase, power, amplitude, normalzied amplitude, db and square amplitude are not affected by the norma vairable.

If Two-Sided (2) is selected for Spectrum Type (st) and Normalize (norma) is set to true, the complex, real, imaginary, magnitude and square magnitude results will be divided by n\,, where n\, is the size of the input signal.

If One-Sided (1) is selected for Spectrum Type (st) and Normalize (norma) is set to true, the complex, real, imaginary, magnitude and square magnitude results will be normalized as follows. Let res_s'\, be the normalized result:

res_s'(i) = 
\begin{cases} 
res_s(i)/n,  & \mbox{if } i=0 \\
2*res_s(i)/n,  & \mbox{otherwise} 
\end{cases}

Automatic Computation of Sampling Interval

The automatically computed sampling interval is the average increment of the time sequence, which is usually from the X column associated with the input signal. If there is no associated X column, the row numbers will be used. Note that if Origin fails to get the average increment, the sampling interval will be set to 1.

Frequency

The frequency column is obtained from the sampling intervalFft1 help English files image014.gifand the number of input data points N. The nth frequency datum is given by:

f_n=\frac n{N\Delta t}

If there are N input data points, the frequency domain will also have N points with the maximum frequency, f_{\max }\,\!equal to\frac 1{\Delta t}(1-\frac 1N).If shift result option is not selected, the transform will display from 0 to f_{\max }\,\!. Otherwise, the shifted transform displays from -\frac{f_{\max }}2 to \frac{f_{\max }}2.

Windows

Windows are used to suppress leakage. Different types of windows are defined as follows in Origin.

In the equations below, n is the index of data and N is the total number of the dataset.

Rectangular Window:

w[n]=1\,\!

Welch Window:

w[n]=1-\left( \frac{n-\frac 12(N-1)}{\frac 12(N+1)}\right) ^2

Triangular Window:

odd: w(n)=\frac 2{N+1}(\frac {N+1}2-|n+1-\frac {N+1}2|)
even: w(n)=\frac 2N(\frac N2-|n+1-\frac {N+1}2|)

Bartlett Window:

w(n)=\frac 2{N-1}(\frac{N-1}2-|n-\frac{N-1}2|)

Hanning Window:

w[n]=\frac 12[1-\cos (\frac{2\pi n}{N-1})]

Hamming Window:

w[n]=0.54-0.46\cos (\frac{2\pi n}{N-1})

Blackman Window:

w[n]=0.42-0.5\cos (\frac{2\pi n}{N-1})+0.08\cos (\frac{4\pi n}{N-1})

Gaussian Window:

w[n]=exp(-0.5(Alpha( \frac{2n}{N-1}-1 ))^2) \,\!

Kaiser Window:

w[n]=I(beta*\sqrt{1-(\frac{2n}{N-1}-1)^2}) / I(beta) \,\!