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_real_3d (c06py)

Purpose

nag_sum_fft_real_3d (c06py) computes the three-dimensional discrete Fourier transform of a trivariate sequence of real data values.

Syntax

[y, ifail] = c06py(n1, n2, n3, x)
[y, ifail] = nag_sum_fft_real_3d(n1, n2, n3, x)

Description

nag_sum_fft_real_3d (c06py) computes the three-dimensional discrete Fourier transform of a trivariate sequence of real data values x j1 j2 j3 x j1 j2 j3 , for j1 = 0,1,,n11j1=0,1,,n1-1, j2 = 0,1,,n21j2=0,1,,n2-1 and j3 = 0,1,,n31j3=0,1,,n3-1.
The discrete Fourier transform is here defined by
n11
k1 k2 k3 = 1/( sqrt( n1 n2 n3 ) ) j2 = 0n21 j3 = 0n31 x j1 j2 j3 × exp(2 π i(( j1 k1 )/(n1) + ( j2 k2 )/(n2) + ( j3 k3 )/(n3))),
j1 = 0
z^ k1 k2 k3 = 1 n1 n2 n3 j1=0 n1-1 j2=0 n2-1 j3=0 n3-1 x j1 j2 j3 × exp( -2πi ( j1 k1 n1 + j2 k2 n2 + j3 k3 n3 ) ) ,
where k1 = 0,1,, n11 k1 = 0,1,, n1-1 , k2 = 0,1,, n21 k2 = 0,1,, n2-1  and k3 = 0,1,, n31 k3 = 0,1,, n3-1 . (Note the scale factor of 1/(sqrt( n1 n2 n3 )) 1 n1 n2 n3  in this definition.)
The transformed values k1 k2 k3 z^ k1 k2 k3  are complex. Because of conjugate symmetry (i.e., k1 k2 k3 z^ k1 k2 k3  is the complex conjugate of (n1k1) k2 k3 z^ ( n1 - k1 ) k2 k3 ), only slightly more than half of the Fourier coefficients need to be stored in the output.
A call of nag_sum_fft_real_3d (c06py) followed by a call of nag_sum_fft_hermitian_3d (c06pz) will restore the original data.
This function performs multiple one-dimensional discrete Fourier transforms by the fast Fourier transform (FFT) algorithm in Brigham (1974) and Temperton (1983).

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:     n1 – int64int32nag_int scalar
n1n1, the first dimension of the transform.
Constraint: n11n11.
2:     n2 – int64int32nag_int scalar
n2n2, the second dimension of the transform.
Constraint: n21n21.
3:     n3 – int64int32nag_int scalar
n3n3, the third dimension of the transform.
Constraint: n31n31.
4:     x( n1 × n2 × n3 n1×n2×n3 ) – double array
The real input dataset xx, where x j1 j2 j3 x j1 j2 j3  is stored in x( j3 × n1 n2 + j2 × n1 + j1 + 1 ) x j3 × n1 n2 + j2 × n1 + j1 +1 , for j1 = 0,1,,n11j1=0,1,,n1-1, j2 = 0,1,,n21j2=0,1,,n2-1 and j3 = 0,1,,n31j3=0,1,,n3-1. That is, if x is regarded as a three-dimensional array of dimension (0 : n11,0 : n21,0 : n31) (0:n1-1,0:n2-1,0:n3-1) , then x(j1,j2,j3)xj1j2j3 must contain xj1j2j3xj1j2j3.

Optional Input Parameters

None.

Input Parameters Omitted from the MATLAB Interface

None.

Output Parameters

1:     y( (n1 / 2 + 1) × n2 × n3 (n1/2+1)×n2×n3 ) – complex array
The complex output dataset z^, where k1 k2 k3 z^ k1 k2 k3 is stored in y( k3 × (n1 / 2 + 1) n2 + k2 × (n1 / 2 + 1) + k1 + 1 ) y k3 × ( n1 /2+1 ) n2 + k2 × ( n1 /2+1) + k1 +1 , for k1 = 0,1,,n1 / 2k1=0,1,,n1/2, k2 = 0,1,,n21k2=0,1,,n2-1 and k3 = 0,1,,n31k3=0,1,,n3-1. That is, if y is regarded as a three-dimensional array of dimension (0 : n1 / 2,0 : n21,0 : n31) (0:n1/2,0:n2-1,0:n3-1) , then y(k1,k2,k3)yk1k2k3 contains k1k2k3z^k1k2k3. Note the first dimension is cut roughly by half to remove the redundant information due to conjugate symmetry.
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
Constraint: n11n11.
  ifail = 2ifail=2
Constraint: n21n21.
  ifail = 3ifail=3
Constraint: n31n31.
  ifail = 4ifail=4
An internal error has occurred in this function. Check the function call and any array sizes. If the call is correct then please contact NAG for assistance.
  ifail = 999ifail=-999
Dynamic memory allocation failed.

Accuracy

Some indication of accuracy can be obtained by performing a forward transform using nag_sum_fft_real_3d (c06py) and a backward transform using nag_sum_fft_hermitian_3d (c06pz), and comparing the results with the original sequence (in exact arithmetic they would be identical).

Further Comments

The time taken by nag_sum_fft_real_3d (c06py) is approximately proportional to n1 n2 n3 log(n1n2n3) n1 n2 n3 log( n1 n2 n3 ) , but also depends on the factors of n1n1, n2n2 and n3n3. nag_sum_fft_real_3d (c06py) is fastest if the only prime factors of n1n1, n2n2 and n3n3 are 22, 33 and 55, and is particularly slow if one of the dimensions is a large prime, or has large prime factors.
Workspace is internally allocated by nag_sum_fft_real_3d (c06py). The total size of these arrays is approximately proportional to n1 n2 n3 n1 n2 n3 .

Example

function nag_sum_fft_real_3d_example
n1 = int64(3);
n2 = int64(3);
n3 = int64(4);
x = zeros(n1, n2, n3);
x(1, :, :) = [1.541, 0.161, 1.989, 0.037;
              0.346, 1.907, 0.001, 1.915;
              1.754, 0.042, 1.991, 0.151];

x(2, :, :) = [0.584, 1.004, 1.408, 0.252;
              1.284, 1.137, 0.467, 1.834;
              0.855, 0.725, 1.647, 0.096];

x(3, :, :) = [0.010, 1.844, 0.452, 1.154;
              1.960, 0.240, 1.424, 0.987;
              0.089, 1.660, 0.708, 0.872];
% Compute Transform.  Reshape x into a 1-d array
[y, ifail] = nag_sum_fft_real_3d(n1, n2, n3, reshape(x,n1*n2*n3,1));

% Display y as a 3-d array
fprintf('\nComponents of discrete Fourier transform:\n');
reshape(y,n1/2,n2,n3)

% Compute Inverse Transform
[x1, ifail] = nag_sum_fft_hermitian_3d(n1, n2, n3, y) ;
fprintf('Original sequence as restored by inverse transform:\n');
reshape(x1, n1, n2, n3)
 

Components of discrete Fourier transform:

ans(:,:,1) =

   5.7547 + 0.0000i  -0.2683 - 0.4203i  -0.2683 + 0.4203i
   0.0814 + 0.0154i   0.0382 + 0.1983i   0.0674 - 0.1220i


ans(:,:,2) =

  -0.2773 - 0.2370i   0.1085 - 0.7560i  -0.6882 + 0.2100i
   0.0605 + 0.1563i  -0.2752 + 0.2954i   0.2804 + 0.0122i


ans(:,:,3) =

   0.4153 + 0.0000i   0.1753 + 0.8712i   0.1753 - 0.8712i
   0.6446 - 0.4779i   1.5848 + 0.6163i  -0.1134 - 1.5552i


ans(:,:,4) =

  -0.2773 + 0.2370i  -0.6882 - 0.2100i   0.1085 + 0.7560i
   0.0469 - 0.0772i   0.2005 + 0.0614i  -0.1281 - 0.1173i

Original sequence as restored by inverse transform:

ans(:,:,1) =

    1.5410    0.3460    1.7540
    0.5840    1.2840    0.8550
    0.0100    1.9600    0.0890


ans(:,:,2) =

    0.1610    1.9070    0.0420
    1.0040    1.1370    0.7250
    1.8440    0.2400    1.6600


ans(:,:,3) =

    1.9890    0.0010    1.9910
    1.4080    0.4670    1.6470
    0.4520    1.4240    0.7080


ans(:,:,4) =

    0.0370    1.9150    0.1510
    0.2520    1.8340    0.0960
    1.1540    0.9870    0.8720


function c06py_example
n1 = int64(3);
n2 = int64(3);
n3 = int64(4);
x = zeros(n1, n2, n3);
x(1, :, :) = [1.541, 0.161, 1.989, 0.037;
              0.346, 1.907, 0.001, 1.915;
              1.754, 0.042, 1.991, 0.151];

x(2, :, :) = [0.584, 1.004, 1.408, 0.252;
              1.284, 1.137, 0.467, 1.834;
              0.855, 0.725, 1.647, 0.096];

x(3, :, :) = [0.010, 1.844, 0.452, 1.154;
              1.960, 0.240, 1.424, 0.987;
              0.089, 1.660, 0.708, 0.872];
% Compute Transform.  Reshape x into a 1-d array
[y, ifail] = c06py(n1, n2, n3, reshape(x,n1*n2*n3,1));

% Display y as a 3-d array
fprintf('\nComponents of discrete Fourier transform:\n');
reshape(y,n1/2,n2,n3)

% Compute Inverse Transform
[x1, ifail] = c06pz(n1, n2, n3, y) ;
fprintf('Original sequence as restored by inverse transform:\n');
reshape(x1, n1, n2, n3)
 

Components of discrete Fourier transform:

ans(:,:,1) =

   5.7547 + 0.0000i  -0.2683 - 0.4203i  -0.2683 + 0.4203i
   0.0814 + 0.0154i   0.0382 + 0.1983i   0.0674 - 0.1220i


ans(:,:,2) =

  -0.2773 - 0.2370i   0.1085 - 0.7560i  -0.6882 + 0.2100i
   0.0605 + 0.1563i  -0.2752 + 0.2954i   0.2804 + 0.0122i


ans(:,:,3) =

   0.4153 + 0.0000i   0.1753 + 0.8712i   0.1753 - 0.8712i
   0.6446 - 0.4779i   1.5848 + 0.6163i  -0.1134 - 1.5552i


ans(:,:,4) =

  -0.2773 + 0.2370i  -0.6882 - 0.2100i   0.1085 + 0.7560i
   0.0469 - 0.0772i   0.2005 + 0.0614i  -0.1281 - 0.1173i

Original sequence as restored by inverse transform:

ans(:,:,1) =

    1.5410    0.3460    1.7540
    0.5840    1.2840    0.8550
    0.0100    1.9600    0.0890


ans(:,:,2) =

    0.1610    1.9070    0.0420
    1.0040    1.1370    0.7250
    1.8440    0.2400    1.6600


ans(:,:,3) =

    1.9890    0.0010    1.9910
    1.4080    0.4670    1.6470
    0.4520    1.4240    0.7080


ans(:,:,4) =

    0.0370    1.9150    0.1510
    0.2520    1.8340    0.0960
    1.1540    0.9870    0.8720



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