hide long namesshow long names
hide short namesshow short names
Integer type:  int32  int64  nag_int  show int32  show int32  show int64  show int64  show nag_int  show nag_int

PDF version (NAG web site, 64-bit version, 64-bit version)
Chapter Contents
Chapter Introduction
NAG Toolbox

NAG Toolbox: nag_sum_fft_realherm_1d_multi_col (c06pq)

Purpose

nag_sum_fft_realherm_1d_multi_col (c06pq) computes the discrete Fourier transforms of mm sequences, each containing nn real data values or a Hermitian complex sequence stored column-wise in a complex storage format.

Syntax

[x, ifail] = c06pq(direct, n, m, x)
[x, ifail] = nag_sum_fft_realherm_1d_multi_col(direct, n, m, x)

Description

Given mm sequences of nn real data values xjp xjp , for j = 0,1,,n1j=0,1,,n-1 and p = 1,2,,mp= 1,2,,m, nag_sum_fft_realherm_1d_multi_col (c06pq) simultaneously calculates the Fourier transforms of all the sequences defined by
n1
kp = 1/(sqrt(n))xjp × exp(i(2πjk)/n),  k0,1,,n1​ and ​p = 1,2,,m.
j = 0
z^kp = 1n j=0 n-1 xjp × exp( -i 2πjk n ) ,   k0,1,,n-1 ​ and ​ p=1,2,,m .
The transformed values kp z^ k p  are complex, but for each value of pp the kp z^ k p  form a Hermitian sequence (i.e., nkp z^ n-k p  is the complex conjugate of kp z^ k p ), so they are completely determined by mn mn  real numbers (since 0p z^0p  is real, as is n / 2p z^ n/2 p  for nn even).
Alternatively, given mm Hermitian sequences of nn complex data values zjp zjp , this function simultaneously calculates their inverse (backward) discrete Fourier transforms defined by
n1
kp = 1/(sqrt(n))zjp × exp(i(2πjk)/n),  k = 0,1,,n1​ and ​p = 1,2,,m.
j = 0
x^kp = 1n j=0 n-1 zjp × exp( i 2πjk n ) ,   k=0,1,,n-1 ​ and ​ p=1,2,,m .
The transformed values kp x^kp  are real.
(Note the scale factor 1/(sqrt(n)) 1n  in the above definition.)
A call of nag_sum_fft_realherm_1d_multi_col (c06pq) with direct = 'F'direct='F' followed by a call with direct = 'B'direct='B' will restore the original data.
The function uses a variant of the fast Fourier transform (FFT) algorithm (see Brigham (1974)) known as the Stockham self-sorting algorithm, which is described in Temperton (1983). Special coding is provided for the factors 22, 33, 44 and 55.

References

Brigham E O (1974) The Fast Fourier Transform Prentice–Hall
Temperton C (1983) Fast mixed-radix real Fourier transforms J. Comput. Phys. 52 340–350

Parameters

Compulsory Input Parameters

1:     direct – string (length ≥ 1)
If the forward transform as defined in Section [Description] is to be computed, then direct must be set equal to 'F'.
If the backward transform is to be computed then direct must be set equal to 'B'.
Constraint: direct = 'F'direct='F' or 'B''B'.
2:     n – int64int32nag_int scalar
nn, the number of real or complex values in each sequence.
Constraint: n1n1.
3:     m – int64int32nag_int scalar
mm, the number of sequences to be transformed.
Constraint: m1m1.
4:     x( (n + 2) × m (n+2)×m ) – double array
The data must be stored in x as if in a two-dimensional array of dimension (0 : n + 1,1 : m)(0:n+1,1:m); each of the mm sequences is stored in a column of the array. In other words, if the data values of the ppth sequence to be transformed are denoted by xjpxjp, for j = 0,1,,n1j=0,1,,n-1, then:
  • if direct = 'F'direct='F', x((p1) × (n + 2) + j)x((p-1)×(n+2)+j) must contain xjpxjp, for j = 0,1,,n1j=0,1,,n-1 and p = 1,2,,mp=1,2,,m;
  • if direct = 'B'direct='B', x ((p1) × (n + 2) + 2 × k)x ( (p-1) × (n+2) + 2×k ) and x ((p1) × (n + 2) + 2 × k + 1)x ( (p-1) × (n+2) + 2×k+1 ) must contain the real and imaginary parts respectively of kpz^kp, for k = 0,1,,n / 2k=0,1,,n/2 and p = 1,2,,mp=1,2,,m. (Note that for the sequence kpz^kp to be Hermitian, the imaginary part of 0pz^0p, and of n / 2p z^ n/2 p  for nn even, must be zero.)

Optional Input Parameters

None.

Input Parameters Omitted from the MATLAB Interface

work

Output Parameters

1:     x( (n + 2) × m (n+2)×m ) – double array
  • if direct = 'F'direct='F' and x is declared with bounds (0 : n + 1,1 : m)(0:n+1,1:m) then x(2 × k,p)x2×kp and x(2 × k + 1,p)x2×k+1p will contain the real and imaginary parts respectively of kp z^ k p , for k = 0,1,,n / 2k=0,1,,n/2 and p = 1,2,,mp=1,2,,m;
  • if direct = 'B'direct='B' and x is declared with bounds (0 : n + 1,1 : m)(0:n+1,1:m) then x(j,p)xjp will contain xjpxjp, for j = 0,1,,n1j=0,1,,n-1 and p = 1,2,,mp=1,2,,m.
2:     ifail – int64int32nag_int scalar
ifail = 0ifail=0 unless the function detects an error (see [Error Indicators and Warnings]).

Error Indicators and Warnings

Errors or warnings detected by the function:
  ifail = 1ifail=1
On entry,m < 1m<1.
  ifail = 2ifail=2
On entry,n < 1n<1.
  ifail = 3ifail=3
On entry,direct'F'direct'F' or 'B''B'.
  ifail = 4ifail=4
An unexpected error has occurred in an internal call. Check all function calls and array dimensions. Seek expert help.

Accuracy

Some indication of accuracy can be obtained by performing a subsequent inverse transform and comparing the results with the original sequence (in exact arithmetic they would be identical).

Further Comments

The time taken by nag_sum_fft_realherm_1d_multi_col (c06pq) is approximately proportional to nm log(n)nm log(n), but also depends on the factors of nn. nag_sum_fft_realherm_1d_multi_col (c06pq) is fastest if the only prime factors of nn are 22, 33 and 55, and is particularly slow if nn is a large prime, or has large prime factors.

Example

function nag_sum_fft_realherm_1d_multi_col_example
direct = 'F';
n = int64(6);
m = int64(3);
x = [0.3854;
     0.6772;
     0.1138;
     0.6751;
     0.6362;
     0.1424;
     0;
     0;
     0.5417;
     0.2983;
     0.1181;
     0.7255;
     0.8638;
     0.8723;
     0;
     0;
     0.9172;
     0.0644;
     0.6037;
     0.643;
     0.0428;
     0.4815;
     0;
     0];
[xOut, ifail] = nag_sum_fft_realherm_1d_multi_col(direct, n, m, x)
 

xOut =

    1.0737
         0
   -0.1041
   -0.0044
    0.1126
   -0.3738
   -0.1467
         0
    1.3961
         0
   -0.0365
    0.4666
    0.0780
   -0.0607
   -0.1521
         0
    1.1237
         0
    0.0914
   -0.0508
    0.3936
    0.3458
    0.1530
         0


ifail =

                    0


function c06pq_example
direct = 'F';
n = int64(6);
m = int64(3);
x = [0.3854;
     0.6772;
     0.1138;
     0.6751;
     0.6362;
     0.1424;
     0;
     0;
     0.5417;
     0.2983;
     0.1181;
     0.7255;
     0.8638;
     0.8723;
     0;
     0;
     0.9172;
     0.0644;
     0.6037;
     0.643;
     0.0428;
     0.4815;
     0;
     0];
[xOut, ifail] = c06pq(direct, n, m, x)
 

xOut =

    1.0737
         0
   -0.1041
   -0.0044
    0.1126
   -0.3738
   -0.1467
         0
    1.3961
         0
   -0.0365
    0.4666
    0.0780
   -0.0607
   -0.1521
         0
    1.1237
         0
    0.0914
   -0.0508
    0.3936
    0.3458
    0.1530
         0


ifail =

                    0



PDF version (NAG web site, 64-bit version, 64-bit version)
Chapter Contents
Chapter Introduction
NAG Toolbox

© The Numerical Algorithms Group Ltd, Oxford, UK. 2009–2013