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_withdraw_fft_real_1d_multi_rfmt (c06fp)

## Purpose

nag_sum_fft_real_1d_multi_rfmt (c06fp) computes the discrete Fourier transforms of $m$ sequences, each containing $n$ real data values. This function is designed to be particularly efficient on vector processors.

## Syntax

[x, trig, ifail] = c06fp(m, n, x, init, trig)
[x, trig, ifail] = nag_sum_withdraw_fft_real_1d_multi_rfmt(m, n, x, init, trig)

## Description

Given $m$ sequences of $n$ real data values ${x}_{\mathit{j}}^{\mathit{p}}$, for $\mathit{j}=0,1,\dots ,n-1$ and $\mathit{p}=1,2,\dots ,m$, nag_sum_fft_real_1d_multi_rfmt (c06fp) simultaneously calculates the Fourier transforms of all the sequences defined by
 $z^ k p = 1n ∑ j=0 n-1 xjp × exp -i 2πjkn , k= 0, 1, …, n-1 ​ and ​ p= 1,2,…,m .$
(Note the scale factor $\frac{1}{\sqrt{n}}$ in this definition.)
The transformed values ${\stackrel{^}{z}}_{k}^{p}$ are complex, but for each value of $p$ the ${\stackrel{^}{z}}_{k}^{p}$ form a Hermitian sequence (i.e., ${\stackrel{^}{z}}_{n-k}^{p}$ is the complex conjugate of ${\stackrel{^}{z}}_{k}^{p}$), so they are completely determined by $mn$ real numbers (see also the C06 Chapter Introduction).
The discrete Fourier transform is sometimes defined using a positive sign in the exponential term:
 $z^kp = 1n ∑ j=0 n-1 xjp × exp +i 2πjkn .$
To compute this form, this function should be followed by forming the complex conjugates of the ${\stackrel{^}{z}}_{k}^{p}$; that is $x\left(\mathit{k}\right)=-x\left(\mathit{k}\right)$, for $\mathit{k}=\left(n/2+1\right)×m+1,\dots ,m×n$.
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$, $3$, $4$, $5$ and $6$. This function is designed to be particularly efficient on vector processors, and it becomes especially fast as $m$, the number of transforms to be computed in parallel, increases.

## 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:     $\mathrm{m}$int64int32nag_int scalar
$m$, the number of sequences to be transformed.
Constraint: ${\mathbf{m}}\ge 1$.
2:     $\mathrm{n}$int64int32nag_int scalar
$n$, the number of real values in each sequence.
Constraint: ${\mathbf{n}}\ge 1$.
3:     $\mathrm{x}\left({\mathbf{m}}×{\mathbf{n}}\right)$ – double array
The data must be stored in x as if in a two-dimensional array of dimension $\left(1:{\mathbf{m}},0:{\mathbf{n}}-1\right)$; each of the $m$ sequences is stored in a row of the array. In other words, if the data values of the $p$th sequence to be transformed are denoted by ${x}_{j}^{p}$, for $\mathit{j}=0,1,\dots ,n-1$, then the $mn$ elements of the array x must contain the values
 $x01 , x02 ,…, x0m , x11 , x12 ,…, x1m ,…, x n-1 1 , x n-1 2 ,…, x n-1 m .$
4:     $\mathrm{init}$ – string (length ≥ 1)
Indicates whether trigonometric coefficients are to be calculated.
${\mathbf{init}}=\text{'I'}$
Calculate the required trigonometric coefficients for the given value of $n$, and store in the array trig.
${\mathbf{init}}=\text{'S'}$ or $\text{'R'}$
The required trigonometric coefficients are assumed to have been calculated and stored in the array trig in a prior call to one of nag_sum_fft_real_1d_multi_rfmt (c06fp), nag_sum_fft_hermitian_1d_multi_rfmt (c06fq) or nag_sum_fft_complex_1d_multi_rfmt (c06fr). The function performs a simple check that the current value of $n$ is consistent with the values stored in trig.
Constraint: ${\mathbf{init}}=\text{'I'}$, $\text{'S'}$ or $\text{'R'}$.
5:     $\mathrm{trig}\left(2×{\mathbf{n}}\right)$ – double array
If ${\mathbf{init}}=\text{'S'}$ or $\text{'R'}$, trig must contain the required trigonometric coefficients that have been previously calculated. Otherwise trig need not be set.

None.

### Output Parameters

1:     $\mathrm{x}\left({\mathbf{m}}×{\mathbf{n}}\right)$ – double array
The $m$ discrete Fourier transforms stored as if in a two-dimensional array of dimension $\left(1:{\mathbf{m}},0:{\mathbf{n}}-1\right)$. Each of the $m$ transforms is stored in a row of the array in Hermitian form, overwriting the corresponding original sequence. If the $n$ components of the discrete Fourier transform ${\stackrel{^}{z}}_{k}^{p}$ are written as ${a}_{k}^{p}+i{b}_{k}^{p}$, then for $0\le k\le n/2$, ${a}_{k}^{p}$ is contained in ${\mathbf{x}}\left(p,k\right)$, and for $1\le k\le \left(n-1\right)/2$, ${b}_{k}^{p}$ is contained in ${\mathbf{x}}\left(p,n-k\right)$. (See also Real transforms in the C06 Chapter Introduction.)
2:     $\mathrm{trig}\left(2×{\mathbf{n}}\right)$ – double array
Contains the required coefficients (computed by the function if ${\mathbf{init}}=\text{'I'}$).
3:     $\mathrm{ifail}$int64int32nag_int scalar
${\mathbf{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:
${\mathbf{ifail}}=1$
 On entry, ${\mathbf{m}}<1$.
${\mathbf{ifail}}=2$
 On entry, ${\mathbf{n}}<1$.
${\mathbf{ifail}}=3$
 On entry, ${\mathbf{init}}\ne \text{'I'}$, $\text{'S'}$ or $\text{'R'}$.
${\mathbf{ifail}}=4$
Not used at this Mark.
${\mathbf{ifail}}=5$
 On entry, ${\mathbf{init}}=\text{'S'}$ or $\text{'R'}$, but the array trig and the current value of n are inconsistent.
${\mathbf{ifail}}=6$
An unexpected error has occurred in an internal call. Check all function calls and array dimensions. Seek expert help.
${\mathbf{ifail}}=-99$
${\mathbf{ifail}}=-399$
Your licence key may have expired or may not have been installed correctly.
${\mathbf{ifail}}=-999$
Dynamic memory allocation failed.

## 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_real_1d_multi_rfmt (c06fp) is approximately proportional to $nm\mathrm{log}\left(n\right)$, but also depends on the factors of $n$. nag_sum_fft_real_1d_multi_rfmt (c06fp) is fastest if the only prime factors of $n$ are $2$, $3$ and $5$, and is particularly slow if $n$ is a large prime, or has large prime factors.

## Example

This example reads in sequences of real data values and prints their discrete Fourier transforms (as computed by nag_sum_fft_real_1d_multi_rfmt (c06fp)). The Fourier transforms are expanded into full complex form using and printed. Inverse transforms are then calculated by conjugating and calling nag_sum_fft_hermitian_1d_multi_rfmt (c06fq) showing that the original sequences are restored.
```function c06fp_example

fprintf('c06fp example results\n\n');

% 3 real sequences stored as rows
m = int64(3);
n = int64(6);
x = [0.3854  0.6772  0.1138  0.6751  0.6362  0.1424;
0.5417  0.2983  0.1181  0.7255  0.8638  0.8723;
0.9172  0.0644  0.6037  0.6430  0.0428  0.4815];

% Transform to get Hermitian sequences
init = 'Initial';
trig = zeros(2*n,1);
[xt, trig, ifail] = c06fp(m, n, x, init, trig);
disp('Discrete Fourier transforms in Hermitian format:');
disp(xt);

for j = 1:m
zt(j,:) = nag_herm2complex(xt(j,:));
end
title = 'Discrete Fourier transforms in full complex format:';
[ifail] = x04da('General','Non-unit', zt, title);

% Restore data by conjugation and back transform
init = 'Subsequent';
nd = double(n);
xt(1:m,floor(nd/2)+2:n) = -xt(1:m,floor(nd/2)+2:n);
[xr, trig, ifail] = c06fq(m, n, xt, init, trig);

fprintf('\n');
disp('Original data as restored by inverse transform:');
disp(xr);

function [z] = nag_herm2complex(x);
n = size(x,2);
z(1) = complex(x(1));
for j = 2:floor((n-1)/2) + 1
z(j) = x(j) + i*x(n-j+2);
z(n-j+2) = x(j) - i*x(n-j+2);
end
if (mod(n,2)==0)
z(n/2+1) = complex(x(n/2+1));
end
```
```c06fp example results

Discrete Fourier transforms in Hermitian format:
1.0737   -0.1041    0.1126   -0.1467   -0.3738   -0.0044
1.3961   -0.0365    0.0780   -0.1521   -0.0607    0.4666
1.1237    0.0914    0.3936    0.1530    0.3458   -0.0508

Discrete Fourier transforms in full complex format:
1          2          3          4          5          6
1      1.0737    -0.1041     0.1126    -0.1467     0.1126    -0.1041
0.0000    -0.0044    -0.3738     0.0000     0.3738     0.0044

2      1.3961    -0.0365     0.0780    -0.1521     0.0780    -0.0365
0.0000     0.4666    -0.0607     0.0000     0.0607    -0.4666

3      1.1237     0.0914     0.3936     0.1530     0.3936     0.0914
0.0000    -0.0508     0.3458     0.0000    -0.3458     0.0508

Original data as restored by inverse transform:
0.3854    0.6772    0.1138    0.6751    0.6362    0.1424
0.5417    0.2983    0.1181    0.7255    0.8638    0.8723
0.9172    0.0644    0.6037    0.6430    0.0428    0.4815

```