Integer type:  int32  int64  nag_int  show int32  show int32  show int64  show int64  show nag_int  show nag_int

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}_{{j}_{1}{j}_{2}{j}_{3}}$, for j1 = 0,1,,n11$\mathit{j1}=0,1,\dots ,{n}_{1}-1$, j2 = 0,1,,n21$\mathit{j2}=0,1,\dots ,{n}_{2}-1$ and j3 = 0,1,,n31$\mathit{j3}=0,1,\dots ,{n}_{3}-1$.
The discrete Fourier transform is here defined by
 n1 − 1 ẑ k1 k2 k3 = 1/( sqrt( n1 n2 n3 ) ) ∑ ∑ j2 = 0n2 − 1 ∑ j3 = 0n3 − 1 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 ${k}_{1}=0,1,\dots ,{n}_{1}-1$, k2 = 0,1,, n21 ${k}_{2}=0,1,\dots ,{n}_{2}-1$ and k3 = 0,1,, n31 ${k}_{3}=0,1,\dots ,{n}_{3}-1$. (Note the scale factor of 1/(sqrt( n1 n2 n3 )) $\frac{1}{\sqrt{{n}_{1}{n}_{2}{n}_{3}}}$ in this definition.)
The transformed values k1 k2 k3 ${\stackrel{^}{z}}_{{k}_{1}{k}_{2}{k}_{3}}$ are complex. Because of conjugate symmetry (i.e., k1 k2 k3 ${\stackrel{^}{z}}_{{k}_{1}{k}_{2}{k}_{3}}$ is the complex conjugate of (n1k1) k2 k3 ${\stackrel{^}{z}}_{\left({n}_{1}-{k}_{1}\right){k}_{2}{k}_{3}}$), 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
n1${n}_{1}$, the first dimension of the transform.
Constraint: n11${\mathbf{n1}}\ge 1$.
2:     n2 – int64int32nag_int scalar
n2${n}_{2}$, the second dimension of the transform.
Constraint: n21${\mathbf{n2}}\ge 1$.
3:     n3 – int64int32nag_int scalar
n3${n}_{3}$, the third dimension of the transform.
Constraint: n31${\mathbf{n3}}\ge 1$.
4:     x( n1 × n2 × n3 ${\mathbf{n1}}×{\mathbf{n2}}×{\mathbf{n3}}$) – double array
The real input dataset x$x$, where x j1 j2 j3 ${x}_{\mathit{j1}\mathit{j2}\mathit{j3}}$ is stored in x( j3 × n1 n2 + j2 × n1 + j1 + 1 )${\mathbf{x}}\left(\mathit{j3}×{n}_{1}{n}_{2}+\mathit{j2}×{n}_{1}+\mathit{j1}+1\right)$, for j1 = 0,1,,n11$\mathit{j1}=0,1,\dots ,{n}_{1}-1$, j2 = 0,1,,n21$\mathit{j2}=0,1,\dots ,{n}_{2}-1$ and j3 = 0,1,,n31$\mathit{j3}=0,1,\dots ,{n}_{3}-1$. That is, if x is regarded as a three-dimensional array of dimension (0 : n11,0 : n21,0 : n31) $\left(0:{\mathbf{n1}}-1,0:{\mathbf{n2}}-1,0:{\mathbf{n3}}-1\right)$, then x(j1,j2,j3)${\mathbf{x}}\left({j}_{1},{j}_{2},{j}_{3}\right)$ must contain xj1j2j3${x}_{{j}_{1}{j}_{2}{j}_{3}}$.

None.

None.

Output Parameters

1:     y( (n1 / 2 + 1) × n2 × n3 $\left({\mathbf{n1}}/2+1\right)×{\mathbf{n2}}×{\mathbf{n3}}$) – complex array
The complex output dataset $\stackrel{^}{z}$, where k1 k2 k3 ${\stackrel{^}{z}}_{\mathit{k1}\mathit{k2}\mathit{k3}}$ is stored in y( k3 × (n1 / 2 + 1) n2 + k2 × (n1 / 2 + 1) + k1 + 1 ) ${\mathbf{y}}\left(\mathit{k3}×\left({n}_{1}/2+1\right){n}_{2}+\mathit{k2}×\left({n}_{1}/2+1\right)+\mathit{k1}+1\right)$, for k1 = 0,1,,n1 / 2$\mathit{k1}=0,1,\dots ,{n}_{1}/2$, k2 = 0,1,,n21$\mathit{k2}=0,1,\dots ,{n}_{2}-1$ and k3 = 0,1,,n31$\mathit{k3}=0,1,\dots ,{n}_{3}-1$. That is, if y is regarded as a three-dimensional array of dimension (0 : n1 / 2,0 : n21,0 : n31) $\left(0:{\mathbf{n1}}/2,0:{\mathbf{n2}}-1,0:{\mathbf{n3}}-1\right)$, then y(k1,k2,k3)${\mathbf{y}}\left({k}_{1},{k}_{2},{k}_{3}\right)$ contains k1k2k3${\stackrel{^}{z}}_{{k}_{1}{k}_{2}{k}_{3}}$. Note the first dimension is cut roughly by half to remove the redundant information due to conjugate symmetry.
2:     ifail – int64int32nag_int scalar
${\mathrm{ifail}}={\mathbf{0}}$ unless the function detects an error (see [Error Indicators and Warnings]).

Error Indicators and Warnings

Errors or warnings detected by the function:
ifail = 1${\mathbf{ifail}}=1$
Constraint: n11${\mathbf{n1}}\ge 1$.
ifail = 2${\mathbf{ifail}}=2$
Constraint: n21${\mathbf{n2}}\ge 1$.
ifail = 3${\mathbf{ifail}}=3$
Constraint: n31${\mathbf{n3}}\ge 1$.
ifail = 4${\mathbf{ifail}}=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 = 999${\mathbf{ifail}}=-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).

The time taken by nag_sum_fft_real_3d (c06py) is approximately proportional to n1 n2 n3 log(n1n2n3) ${n}_{1}{n}_{2}{n}_{3}\mathrm{log}\left({n}_{1}{n}_{2}{n}_{3}\right)$, but also depends on the factors of n1${n}_{1}$, n2${n}_{2}$ and n3${n}_{3}$. nag_sum_fft_real_3d (c06py) is fastest if the only prime factors of n1${n}_{1}$, n2${n}_{2}$ and n3${n}_{3}$ are 2$2$, 3$3$ and 5$5$, 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 ${n}_{1}{n}_{2}{n}_{3}$.

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

```