help-gsl
[Top][All Lists]

## Re: [Help-gsl] FFT

 From: Francesco Subject: Re: [Help-gsl] FFT Date: Mon, 06 Jan 2014 20:57:58 +0100 User-agent: Mozilla/5.0 (X11; Linux i686; rv:24.0) Gecko/20100101 Thunderbird/24.2.0

Hi Klaus,

On 06/01/2014 15:30, Klaus Huthmacher wrote:

Dear fellows,


I have another question concerning the FFT of the GSL. I wanted to start low and tried to transform a sine signal.


Since I assume a data set I take into account real values and a length of any magnitude. Thus I end up with Section 16.7: Mixed - radix FFT routines for real data, am I right?
Right.


But to be honest I do not fully grasp between the two methods real_transform and halfcomplex_transform. The first leads to an expected \delta - peak while the second method leads to 'oscillating boundaries'.
I recognize that the documentation is not very clear. Actually the real_transform does a direct FFT transform in place by storing the result in half-complex form. The half_complex transform perform the *inverse* FFT transform by accepting data in half-complex format and returning real data.


The GSL documentation show a clear example where some real temporal data are transformed using real_transform, them the half-complex vector is modified to remove some higher frequencies and the the inverse transform is done with halfcomplex_transform. The net result is the removal of high frequencies from the time resolved signal.


Furthermore I am highly confused that the FFT functions only need the data in the sense of the 'values' of a function and not the corresponding x values.
Well, in the FFT transform the data are implicitly supposed to be taken with uniform sampling. These is a consequence of its mathematic definition given in the GSL FFT chapter. For the direct (forward) transform:

x_j = \sum_{k=0}^{n-1} z_k \exp(-2\pi i j k / n)


The fact that the exponential is proportional to "k" means that the temporal data is given with a uniform sampling. You should compare that to the definition of the fourier transform:

x_j = \int_{t = 0}^{1} z(t) \exp(- 2 \pi i j t / n) dt


I thought to be quiete familier with the fouriertransform by solving equations with paper and pencil but this GSL implementation is curious.
The FFT is the Fourier Transform but of a signal taken with a finite, fixed, sampling step. The Fourier Transform from book assume the signal is defined in an interval of real numbers so the time variable is actually continuous. In the FFT the time is "discrete".


Another source of confusion is the packing convention with the halfcomplex format. The FFT routine are conceived to be efficient and this lead to some strange packing forms.


Another important difference is that the FFT does need a finite set of frequencies and the negative frequencies aliases the higher frequencies. These are consequences of the finite sampling and does not happens in the fourier transform over a continuous variable.
Can someone give me a little more insight?

Your example is quite right. If you modify it in this way:

std::vector<double> data;
unsigned N = 512;
unsigned order = 4;

for (unsigned i = 0; i < N; i++)
{
double x = double(order * i) * 2 * M_PI / double(N);
data.push_back( sin( x ) );
}

The output will be a pure "pulse" (dirac's delta). Here what I obtain:

5.04177e-15
-2.7399e-15
7.42739e-15
-3.21873e-16
7.66043e-15
5.63479e-15
3.9062e-14
-1.17761e-13
-256
-3.24479e-15
-4.43331e-14
4.93448e-16
-2.25068e-14
-4.25877e-15
-1.8501e-14

Essentially everything is zero except the term corresponding to "order".

Francesco