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_realherm_1d_multi_col (c06pq)

## Purpose

nag_sum_fft_realherm_1d_multi_col (c06pq) computes the discrete Fourier transforms of m$m$ sequences, each containing n$n$ 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 m$m$ sequences of n$n$ real data values xjp ${x}_{\mathit{j}}^{\mathit{p}}$, for j = 0,1,,n1$\mathit{j}=0,1,\dots ,n-1$ and p = 1,2,,m$\mathit{p}=1,2,\dots ,m$, nag_sum_fft_realherm_1d_multi_col (c06pq) simultaneously calculates the Fourier transforms of all the sequences defined by
 n − 1 ẑkp = 1/(sqrt(n)) ∑ xjp × exp( − i(2πjk)/n),  k0,1, … ,n − 1​ 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 ${\stackrel{^}{z}}_{k}^{p}$ are complex, but for each value of p$p$ the kp ${\stackrel{^}{z}}_{k}^{p}$ form a Hermitian sequence (i.e., nkp ${\stackrel{^}{z}}_{n-k}^{p}$ is the complex conjugate of kp ${\stackrel{^}{z}}_{k}^{p}$), so they are completely determined by mn $mn$ real numbers (since 0p ${\stackrel{^}{z}}_{0}^{p}$ is real, as is n / 2p ${\stackrel{^}{z}}_{n/2}^{p}$ for n$n$ even).
Alternatively, given m$m$ Hermitian sequences of n$n$ complex data values zjp ${z}_{j}^{p}$, this function simultaneously calculates their inverse (backward) discrete Fourier transforms defined by
 n − 1 x̂kp = 1/(sqrt(n)) ∑ zjp × exp(i(2πjk)/n),  k = 0,1, … ,n − 1​ 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 ${\stackrel{^}{x}}_{k}^{p}$ are real.
(Note the scale factor 1/(sqrt(n)) $\frac{1}{\sqrt{n}}$ in the above definition.)
A call of nag_sum_fft_realherm_1d_multi_col (c06pq) with direct = 'F'${\mathbf{direct}}=\text{'F'}$ followed by a call with direct = 'B'${\mathbf{direct}}=\text{'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 2$2$, 3$3$, 4$4$ and 5$5$.

## 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'${\mathbf{direct}}=\text{'F'}$ or 'B'$\text{'B'}$.
2:     n – int64int32nag_int scalar
n$n$, the number of real or complex values in each sequence.
Constraint: n1${\mathbf{n}}\ge 1$.
3:     m – int64int32nag_int scalar
m$m$, the number of sequences to be transformed.
Constraint: m1${\mathbf{m}}\ge 1$.
4:     x( (n + 2) × m $\left({\mathbf{n}}+2\right)×{\mathbf{m}}$) – double array
The data must be stored in x as if in a two-dimensional array of dimension (0 : n + 1,1 : m)$\left(0:{\mathbf{n}}+1,1:{\mathbf{m}}\right)$; each of the m$m$ sequences is stored in a column of the array. In other words, if the data values of the p$p$th sequence to be transformed are denoted by xjp${x}_{\mathit{j}}^{p}$, for j = 0,1,,n1$\mathit{j}=0,1,\dots ,n-1$, then:
• if direct = 'F'${\mathbf{direct}}=\text{'F'}$, x((p1) × (n + 2) + j)${\mathbf{x}}\left(\left(\mathit{p}-1\right)×\left({\mathbf{n}}+2\right)+\mathit{j}\right)$ must contain xjp${x}_{\mathit{j}}^{\mathit{p}}$, for j = 0,1,,n1$\mathit{j}=0,1,\dots ,n-1$ and p = 1,2,,m$\mathit{p}=1,2,\dots ,m$;
• if direct = 'B'${\mathbf{direct}}=\text{'B'}$, x ((p1) × (n + 2) + 2 × k)${\mathbf{x}}\left(\left(\mathit{p}-1\right)×\left({\mathbf{n}}+2\right)+2×\mathit{k}\right)$ and x ((p1) × (n + 2) + 2 × k + 1)${\mathbf{x}}\left(\left(\mathit{p}-1\right)×\left({\mathbf{n}}+2\right)+2×\mathit{k}+1\right)$ must contain the real and imaginary parts respectively of kp${\stackrel{^}{z}}_{\mathit{k}}^{\mathit{p}}$, for k = 0,1,,n / 2$\mathit{k}=0,1,\dots ,n/2$ and p = 1,2,,m$\mathit{p}=1,2,\dots ,m$. (Note that for the sequence kp${\stackrel{^}{z}}_{k}^{p}$ to be Hermitian, the imaginary part of 0p${\stackrel{^}{z}}_{0}^{p}$, and of n / 2p${\stackrel{^}{z}}_{n/2}^{p}$ for n$n$ even, must be zero.)

None.

work

### Output Parameters

1:     x( (n + 2) × m $\left({\mathbf{n}}+2\right)×{\mathbf{m}}$) – double array
• if direct = 'F'${\mathbf{direct}}=\text{'F'}$ and x is declared with bounds (0 : n + 1,1 : m)$\left(0:{\mathbf{n}}+1,1:{\mathbf{m}}\right)$ then x(2 × k,p)${\mathbf{x}}\left(2×\mathit{k},\mathit{p}\right)$ and x(2 × k + 1,p)${\mathbf{x}}\left(2×\mathit{k}+1,\mathit{p}\right)$ will contain the real and imaginary parts respectively of kp${\stackrel{^}{z}}_{\mathit{k}}^{\mathit{p}}$, for k = 0,1,,n / 2$\mathit{k}=0,1,\dots ,n/2$ and p = 1,2,,m$\mathit{p}=1,2,\dots ,m$;
• if direct = 'B'${\mathbf{direct}}=\text{'B'}$ and x is declared with bounds (0 : n + 1,1 : m)$\left(0:{\mathbf{n}}+1,1:{\mathbf{m}}\right)$ then x(j,p)${\mathbf{x}}\left(\mathit{j},\mathit{p}\right)$ will contain xjp${x}_{\mathit{j}}^{\mathit{p}}$, for j = 0,1,,n1$\mathit{j}=0,1,\dots ,n-1$ and p = 1,2,,m$\mathit{p}=1,2,\dots ,m$.
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$
 On entry, m < 1${\mathbf{m}}<1$.
ifail = 2${\mathbf{ifail}}=2$
 On entry, n < 1${\mathbf{n}}<1$.
ifail = 3${\mathbf{ifail}}=3$
 On entry, direct ≠ 'F'${\mathbf{direct}}\ne \text{'F'}$ or 'B'$\text{'B'}$.
ifail = 4${\mathbf{ifail}}=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).

The time taken by nag_sum_fft_realherm_1d_multi_col (c06pq) is approximately proportional to nm log(n)$nm\mathrm{log}\left(n\right)$, but also depends on the factors of n$n$. nag_sum_fft_realherm_1d_multi_col (c06pq) is fastest if the only prime factors of n$n$ are 2$2$, 3$3$ and 5$5$, and is particularly slow if n$n$ 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

Chapter Contents
Chapter Introduction
NAG Toolbox

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