SFT

SFT is the discrete implementation of the SinFT transform.

At some fixed number \(N\) of nodes

the SFT of array \(f\) of length \(N\!-\!1\) appears as \(g=\,\)SFT\((f)\) defined with

\(\displaystyle g_m=\sum_{n=1}^{N-1} \sin\left( \frac{\pi}{N}\, m\, n\right) \, f_n\)

For \(N\!=\!2^n\) for some integer \(n\), there exist efficient numerical algorithms for evaluation of this sum.

SFT and CFT

SFT and its analogy CFT (for the CosFT transform) are implemented in C. The fast algorithm for evaluation of SFT or CRB performs the evaluation taking of order to \(N \log_2(N)\) operations. One of implementation of the fast algorithms is suggested by the Numerical recipes in C. .

The double version of this algorithm is stored as scft.cin; the name means that SFT and CFT are combined in the single file.

Testing

The simple test of SFT includes encoding of several numbers with amplitudes of the corresponding sinusoidals, and the recovery of the coefficients with the SFT algorithm. File scft.cin should be loaded in order to compile the test below:

#include <stdio.h>
#include <math.h>
#include "scft.cin"
#define NP 16
int main(){ int i;double a[NP],b[NP];
for(i=0;i<NP;i++){ a[i]=
          1.*sin((M_PI/NP)*i))
        + 2.*sin( 2*(M_PI/NP)*i)
        + 3.*sin( 3*(M_PI/NP)*i)
        + 4.*sin( 4*(M_PI/NP)*i)
        + 5.*sin( 5*(M_PI/NP)*i)
        + 6.*sin( 6*(M_PI/NP)*i)
        + 7.*sin( 7*(M_PI/NP)*i)
        + 8.*sin( 8*(M_PI/NP)*i)
        + 9.*sin( 9*(M_PI/NP)*i)
        +10.*sin(10*(M_PI/NP)*i)
        +11.*sin(11*(M_PI/NP)*i)
        +12.*sin(12*(M_PI/NP)*i)
        +13.*sin(13*(M_PI/NP)*i)
        +14.*sin(14*(M_PI/NP)*i)
        +15.*sin(15*(M_PI/NP)*i)
        ;
        b[i]=a[i]/8.;}
sinft(b-1,NP);
for(i=0;i<NP;i++) printf("%2d %19.14lf %19.14lf\n",i,a[i],b[i]);
}
//

The test provides the following output:

 0 0.00000000000000 0.00000000000000
 1 81.22536310087088 1.00000000000000
 2 -40.21871593700679 1.99999999999999
 3 26.37246567150657 2.99999999999999
 4 -19.31370849898477 4.00000000000001
 5 14.96694729431511 5.00000000000000
 6 -11.97284610132397 6.00000000000000
 7 9.74802820470383 6.99999999999999
 8 -7.99999999999999 8.00000000000001
 9 6.56543032662929 9.00000000000000
10 -5.34542910335433 10.00000000000000
11 4.27608908760625 11.00000000000001
12 -3.31370849898468 11.99999999999998
13 2.42677346885870 12.99999999999998
14 -1.59129893903725 14.00000000000004
15 0.78793122685728 14.99999999999998

As it is seen from the example above, the correct call with array "data" of NP points, from zero to NP-1, is

sinft(data-1,NP);

The result is placed at the same array; so, if the original array is also necessary, it may have sense to work not with the original array, but with its copy. This copy may be scaled with factor 2/NP; then, the double call of the same function returns the original array with small rounding errors. For NP=16, these errors are of order of \(10^{-15}\) of the magnitude of numbers in the array.

Characters "-1" appear in the call of the routine. There is deep philosophy and, may be, even religion about representation of arrays in C++ and Fortran, described in the book Numerical recipes in C, with respect to these characters. Bu default, if one omits these two characters in the call of the routine, the significant date begins not from the 0th element of the array, and even not from the 1st element, but from the second one. For the practical applications, one may consider this "-1" as and abracadabra and just not to forget to type it in the programming.

Keywords

C++, CosFT, FFT, Integral transform, SinFT,