Функция генерирует множители таблицы вращения (radix-2 FFT) FFT для функций FFT по основанию 2 (radix-2 FFT).

#include < filter.h >

void twidfftrad2_fr16(complex_fract16 twiddle_table[],
                      int fft_size);

void twidfftrad2_fr32(complex_fract32 twiddle_table[],
                      int fft_size);

Функции twidfftrad2 вычисляют комплексные коэффициенты вращения для radix-2 FFT в таблице размером fft_size, которые будут возвращены в векторе twiddle_table. Размер 

этого вектора, который также называют таблицей вращения (twiddle table), должен быть как минимум fft_size/2. Здесь содержатся пары значений синуса и косинуса, которые используются функцией FFT для вычисления алгоритма БПФ (Fast Fourier Transform, FFT). Таблица, генерируемая функцией twidfftrad2_fr16, может использоваться любой из функций cfft_fr16, ifft_fr16, rfft_fr16 и rfft_fx_fr16, и таблица, генерируемая функцией twidfftrad2_fr32, может использоваться любой из функций cfft_fr32, ifft_fr32, rfft_fr32 и rfft_fx_fr32.

Количество точек FFT fft_size должно быть числом, равным степени 2, и не меньше 8.

Таблица вращения указанного размера содержит значения констант, так что такая таблица обычно генерируется один раз в процессе разработки приложения, после чего она может быть сохранена в памяти в любом удобном, готовом виде (т. е. в реальном приложении будет готовая таблица вращения, без использования функции twidfftrad2).

Приложению, которое вычисляет преобразования FFT разных размеров, не требуется использовать несколько таблиц вращения. Можно использовать единственную таблицу вращения максимального размера (из имеющихся FFT) для вычисления всех размеров FFT. Для этой цели каждая из FFT-функций cfft, ifft и rfft имеют аргумент шага таблицы вращения (twiddle_stride), который приложение должно установить в 1, когда генерируется FFT с самым большим количеством точек данных. Для генерации FFT меньшего размера, этот аргумент шага должен быть установлен по следующей формуле:                      

                  самый большой размер FFT
twiddle_stride = -------------------------
                    текущий размер FFT

Например, если таблица вращения была создана для FFT на 1024 точек, то та же самая таблица вращения может использоваться и для FFT на 256 точек путем установки аргумента шага twiddle_stride = 4.

[Алгоритм]

Функция принимает длину FFT во входном параметре fft_size, и генерирует таблицу комплексных коэффициентов вращения. Вещественная и мнимая части значений в таблице генерируются по формулам:

                    2pi
twid_re(k) =  cos (-----k)
                     n

                    2pi
twid_im(k) = -sin (-----k)
                     n

где

n = fft_size,
k = {0, 1, 2, ..., n/2 - 1}.

[Пример]

#include < filter.h >

#define FFT_SIZE1 256

#define FFT_SIZE2 64

#define TWID_SIZE (FFT_SIZE1/2)

complex_fract32 input1[FFT_SIZE1];
complex_fract32 output1[FFT_SIZE1];
complex_fract32 input2[FFT_SIZE2];
complex_fract32 output2[FFT_SIZE2];
complex_fract32 twiddles[TWID_SIZE];

int block_exponent1, block_exponent2;

int scale_method = 1; 

twidfftrad2_fr32 (twiddles, FFT_SIZE1);
cfft_fr32 (input1,
           output1,
           twiddles,
           1,
           FFT_SIZE1,
           &block_exponent1,
           scale_method);
cfft_fr32 (input1,
           output2,
           twiddles,
           (FFT_SIZE1/FFT_SIZE2),
           FFT_SIZE2,
           &block_exponent2,
           scale_method);

Функция генерирует таблицу множителей вращения для БПФ по основанию 4 (FFT twiddle factors for radix-4 FFT).

#include < filter.h >

void twidfftrad4_fr16(complex_fract16 twiddle_table[],
                      int fft_size);

void twidfft_fr16(complex_fract16 twiddle_table[],
                  int fft_size);

Функция twidfftrad4_fr16 инициализирует таблицу комплексных множителей вращения для БПФ по основанию 4. Количество точек в FFT задается параметром fft_size, и коэффициенты будут возвращены в таблице twiddle_table.

Размер таблицы вращения должен быть как минимум 3*fft_size/4. Таблица может использоваться для нескольких FFT разных размеров путем выделения таблицы максимального размера, и затем использования аргумента шага по таблице (twiddle_stride) для функции FFT. Если шаг равен 1, то функция FFT будет использовать всю таблицу; если FFT имеет размер 1/4 от максимального размера FFT, то тогда шаг должен быть установлен равным 4.

Количество точек FFT fft_size должно быть числом, равным степени 2, и не меньше 16.

Для эффективности таблица вращения обычно генерируется только один раз во время инициализации программы, после чего эта таблица может использоваться в программе любое количество раз как аргумент функции FFT.

Примечание: функция twidfftrad4_fr16 и функции radix-4 FFT предоставлены только в целях совместимости с уже существующими программами. Новые приложения должны использовать вместо этого одну из функций radix-2 FFT (см. cfft, ifft, rfft). Таблица вращения для функций radix-2 FFT может быть сгенерирована вызовом функции twidfftrad2_fr16.

Функция twidfft_fr16 может использоваться как альтернатива функции twidfftrad4_fr16. Они обе работают одинаково.

[Алгоритм]

Функция принимает длину FFT во входном параметре fft_size, и генерирует таблицу комплексных коэффициентов вращения. Вещественная и мнимая части значений в таблице генерируются по формулам:

                    2pi
twid_re(k) =  cos (-----k)
                     n

                    2pi
twid_im(k) = -sin (-----k)
                     n

где

n = fft_size,
k = {0, 1, 2, ..., (3/4)n-1}.

Функция генерирует таблицу множителей вращения для ускоренного алгоритма БПФ.

#include < filter.h >

void twidfftf_fr16(complex_fract16 twiddle_table[],
                   int fft_size);

void twidfftf_fr32(complex_fract32 twiddle_table[],
                   int fft_size);

Функция twidfftf_fr16 генерирует таблицу вращения из комплексных множителей, предназначенную для быстрой версии БПФ по основанию 4 в функции cfftf_fr16, в то время как twidfftf_fr32 генерирует аналогичные множители для функций с ускоренным алгоритмом БПФ смешанного основания cfftf_fr32, ifftf_fr32, rfftf_fr32 и rfftf_fx_fr32. Множители вращения это пары значений косинусов и синусов, сохраняемые в векторе twiddle_table; функции FFT используют эту таблицу для генерации быстрого преобразования Фурье (БПФ, Fast Fourier Transform, FFT). Размер таблицы вращения должен быть как минимум 3*fft_size/4 от fft_size (количества точек в FFT).

Количество точек в FFT должно быть степенью числа 4, и не менее 16 для функции cfftf_fr16, и степенью числа 2 и не меньше 16 для функций cfftf_fr32, ifftf_fr32, rfftf_fr32 и rfftf_fx_fr32.

Таблица вращения указанного размера содержит значения констант, так что такая таблица обычно генерируется один раз в процессе разработки приложения, после чего она может быть сохранена в памяти в любом удобном, готовом виде (т. е. в реальном приложении будет готовая таблица вращения, без использования функции twidfftf).

Приложению, которое вычисляет преобразования FFT разных размеров, не требуется использовать несколько таблиц вращения. Можно использовать единственную таблицу вращения максимального размера (из имеющихся FFT) для вычисления всех размеров FFT. Для этой цели каждая из FFT-функций cfft, ifft и rfft имеют аргумент шага таблицы вращения (twiddle_stride), который приложение должно установить в 1, когда генерируется FFT с самым большим количеством точек данных. Для генерации FFT меньшего размера, этот аргумент шага должен быть установлен по следующей формуле:                      

                  самый большой размер FFT
twiddle_stride = -------------------------
                    текущий размер FFT

Например, если таблица вращения была создана для FFT на 1024 точек, то та же самая таблица вращения может использоваться и для FFT на 256 точек путем установки аргумента шага twiddle_stride = 4.

[Алгоритм]

Функция принимает длину FFT во входном параметре fft_size, и генерирует таблицу комплексных коэффициентов вращения. Вещественная и мнимая части значений в таблице генерируются по формулам:

                    2pi
twid_re(k) =  cos (-----k)
                     n

                    2pi
twid_im(k) = -sin (-----k)
                     n

где

n = fft_size,
k = {0, 1, 2, ..., (3/4)n-1}.

[Пример]

#include < filter.h >

#define FFT_SIZE1 256

#define FFT_SIZE2 64

#define TWIDDLE_SIZE ((3*FFT_SIZE1)/4)

complex_fract32 in1[FFT_SIZE1];
complex_fract32 out1[FFT_SIZE1];
complex_fract32 in2[FFT_SIZE2];
complex_fract32 out2[FFT_SIZE2];
complex_fract32 twiddles[TWIDDLE_SIZE];

 
twidfftf_fr32 (twiddles, FFT_SIZE1);
cfftf_fr32(in1, out1, twiddles, 1, FFT_SIZE1);
cfftf_fr32(in2, out2, twiddles, FFT_SIZE1/FFT_SIZE2, FFT_SIZE2);

Функция для вычисления БПФ по основанию 2 на N точках входного сигнала, представленного вещественными числами (N-point radix-2 real input FFT).

#include < filter.h >

void rfft_fr16(const fract16 input[],
               complex_fract16 output[],
               const complex_fract16 twiddle_table[],
               int twiddle_stride,
               int fft_size,
               int *block_exponent,
               int scale_method);

void rfft_fx_fr16(const _Fract input[],
                  complex_fract16 output[],
                  const complex_fract16 twiddle_table[],
                  int twiddle_stride,
                  int fft_size,
                  int *block_exponent,
                  int scale_method);

void rfft_fr32(const fract32 input[],
               complex_fract32 output[],
               const complex_fract32 twiddle_table[],
               int twiddle_stride,
               int fft_size,
               int *block_exponent,
               int scale_method);

void rfft_fx_fr32(const long _Fract input[],
                  complex_fract32 output[],
                  const complex_fract32 twiddle_table[],
                  int twiddle_stride,
                  int fft_size,
                  int *block_exponent,
                  int scale_method);

Функции rfft выполняют трансформацию последовательности входного сигнала вещественных значений из домена времени в частотный домен с использованием БПФ по основанию 2 (radix-2 FFT). Достоинство этих функций в том, что мнимая часть входного сигнала равна 0, что устраняет половину умножений в алгоритме бабочки.

Размер входного массива (параметр input) и выходного массива (параметр output) равен fft_size, где fft_size представляет количество точек в FFT. Если входные данные могут быть перезаписаны, то оптимального использования памяти можно достичь путем указания входного массива как выходного, при этом размер входного массива должен быть как минимум 2*fft_size.

Таблица вращения передается в аргументе twiddle_table, и в этой таблице должно содержаться как минимум fft_size/2 множителей вращения. Эта таблица составлена из коэффициентов косинуса и синуса (+cosine и -sine), и она может быть инициализирована функцией twidfftrad2_fr16 для использования с функциями rfft_fr16 или rfft_fx_fr16, и twidfftrad2_fr32 для использования с функциями rfft_fr32 или rfft_fx_fr32. Для оптимальной производительности таблица вращения должна быть размещена в другой физической секции памяти, чем выходной массив (параметр output).

Аргумент twiddle_stride (шаг таблицы вращения) должен быть установлен в 1, если таблица вращения была изначально создана для FFT размером fft_size. Если таблица вращения была создана для большего размера FFT (N*fft_size, где N является числом, равным степени 2), то twiddle_stride должен быть установлен в N. Таким образом, этот аргумент предоставляет способ использовать одну таблицу вращения для вычисления преобразований FFT различных размеров.

Аргумент scale_method управляет тем, как какое функция применит масштабирование при вычислении преобразования Фурье. Доступные опции - статическое масштабирование (деление входа на 2 на всех стадиях вычислений), динамическое масштабирование (делением входа на 2 на любой стадии, если самое большое абсолютное входное значение больше или равно 0.25), или без масштабирования. Имейте в виду, что количество стадий, необходимое для вычисления FFT, зависит от размера FFT, и может быть получено по формуле log2(fft_size).

Если выбрано статическое масштабирование, то функция всегда будет выполнять промежуточное масштабирование, предотвращая тем самым переполнение. Потеря точности увеличивается с ростом fft_size, что больше всего влияет на входные сигналы с малой магнитудой (поскольку выходной сигнал будет масштабирован на множитель 1/fft_size). Чтобы выбрать статическое масштабирование, установить аргумент scale_method в значение 1. Будет возвращено значение block_exponent, равное log2(fft_size).

Если выбрано динамическое масштабирование, то функция будет инспектировать промежуточные результаты и применять масштабирования только для случаев, когда возможно переполнение. Потеря точности увеличивается с ростом fft_size, что больше всего влияет на входные сигналы с большой амплитудой (поскольку в этих условиях повышается необходимость в масштабировании). Необходимость инспектирования промежуточных результатов отрицательно влияет на производительность. Чтобы выбрать динамическое масштабирование, установите scale_method в значение 2. Будет возвращено значение block_exponent в диапазоне между 0 и log2(fft_size), в зависимости от того, сколько раз было выполнено масштабирование промежуточных результатов в процессе вычисления.

Если не выбрано масштабирование, то функция никогда не будет масштабировать вниз промежуточные результаты. В этом случае не будет потери точности, если не было переполнения (при переполнении функция вернет результаты с насыщением). Вероятность появления насыщения растет с ростом fft_size, что больше всего актуально при работе с сигналами большой магнитуды. Чтобы выключить масштабирование, установить scale_method в значение 3. Будет возвращено block_exponent, равное 0.

Примечание: любое значение аргумента scale_method, отличающееся от 2 или 3, приведет к выполнению функцией статического масштабирования.

[Условия возникновения ошибки]

Функции rfft немедленно оборвут вычисления, если размер FFT (параметр fft_size) меньше 8, или если параметр twiddle_stride меньше 1.

Длина входной последовательности домена времени fft_size должна быть числом, равным степени 2, и быть не меньше 8.

[Алгоритм]

См. описание функции cfft.

[Пример]

#include < filter.h >

#define FFT_SIZE1 32

#define FFT_SIZE2 256

#define TWID_SIZE (FFT_SIZE2/2)

fract32 in1[FFT_SIZE1], in2[FFT_SIZE2];
complex_fract32 out1[FFT_SIZE1], out2[FFT_SIZE2];
complex_fract32 twiddle[TWID_SIZE];

int block_exponent1, block_exponent2;

twidfftrad2_fr32 (twiddle, FFT_SIZE2);
rfft_fr32 (in1,
           out1,
           twiddle,
           (FFT_SIZE2 / FFT_SIZE1),
           FFT_SIZE1,
           &block_exponent1,
           1 /*статическое масштабирование*/);
rfft_fr32 (in2,
           out2,
           twiddle,
           1,
           FFT_SIZE2,
           &block_exponent2,
           2 /*динамическое масштабирование*/);

Функция для ускоренного вычисления БПФ на N точках входного сигнала, представленного вещественными числами (Fast N-point real input FFT).

#include < filter.h >

void rfftf_fr32(const fract32 input[],
                complex_fract32 output[],
                const complex_fract32 twiddle_table[],
                int twiddle_stride,
                int fft_size);

void rfftf_fx_fr32(const long _Fract input[],
                   complex_fract32 output[],
                   const complex_fract32 twiddle_table[],
                   int twiddle_stride,
                   int fft_size);

Функции rfftf выполняют трансформацию последовательности входного сигнала вещественных значений из домена времени в частотный домен с использованием ускоренной версии БПФ, так называемого быстрого преобразования Фурье (Fast Fourier Transform, FFT). Они выполняют децимацию по частоте, используя алгоритм по смешанному основанию (mixed-radix).

Размер входного массива (параметр input) и выходного массива (параметр output) равен fft_size, где fft_size представляет количество точек в FFT. Количество точек FFT должно быть числом, равным степени двойки, и не меньше 16.

Поскольку комплексный спектр вещественного FFT симметричен относительно средней точки, функции rfftf генерируют только первые (fft_size/2)+1 точек FFT, и размер выходного массива output должен быть как минимум размера (fft_size/2)+1. После завершения функции массив output будет содержать следующие значения:

• Постоянную составляющую сигнала (DC) в ячейке output[0].re (где output[0].im = 0).
• Первую половину комплексного спектра в в ячейках output[1] .. output[(fft_size/2)-1].
• Частоту Найквиста в ячейке output[fft_size/2].re (где output[fft_size/2].im = 0).

Ниже см. секцию "Пример", где показано, как приложение должно конструировать полный комплексный спектр с использование симметрии вещественного FFT.

В аргументе twiddle_table передается ссылка на таблицу вращения, которая должна содержать как минимум 3*fft_size/4 комплексных множителей вращения. Эта таблица должна быть инициализирована комплексными множителями вращения, где вещественная часть коэффициентов содержит положительные значения косинусов, и мнимая часть коэффициентов содержит отрицательные значения синусов. Для инициализации этого массива может использоваться функция twidfftf_fr32.

Если таблица вращения была сгенерирована для fft_size FFT, то аргумент twiddle_stride (шаг таблицы вращения) должен быть установлен в 1. Другими словами, если таблица вращения была сгенерирована для FFT размера x, где x > fft_size, то аргумент twiddle_stride должен быть установлен в x/fft_size. Таким образом, аргумент twiddle_stride позволяет использовать одну и ту же таблицу вращения для разных размеров FFT. При этом аргумент twiddle_stride не может быть нулевым или отрицательным.

Рекомендуется, чтобы выходной массив не был расположен в том же самом подбанке памяти 4K, что и входной массив или таблица вращения, иначе скорость работы функции может снизиться из-за коллизий операций обращения к ячейкам памяти, находящихся в одном банке данных.

Функции используют статическое масштабирование промежуточных результатов, чтобы избежать переполнения, так что выходное значение будет масштабировано вниз на коэффициент 1/fft_size.

[Алгоритм]

Ниже приведена формула базового алгоритма.

x(n) выходной сигнал для отсчета n, где n=0..N-1.

N количество отсчетов дискретного сигнала Y(n).

X(k) амплитуда частотной составляющей, k=0..N-1, т. е. k соответствуют частоты спектра. Частоты = kF/N, где F частота дискретизации.

W_n коэффициенты матрицы вращения.

Набор коэффициентов Xk называется амплитудным спектром сигнала. Частоты синусоид спектра равномерно распределены от 0 (постоянная составляющая) до F/2 (максимально возможная частота в цифровом сигнале).

Реализация функций использует алгоритм смешанного основания (radix-4 и radix-2, т. е. по основанию 4 и по основанию 2).

[Пример]

#include < filter.h >

#include < complex.h >

#define FFTSIZE 32

#define TWIDSIZE ((3 * FFTSIZE) / 4)

fract32 sigdata[FFTSIZE];
complex_fract32 r_output[FFTSIZE];
complex_fract32 twiddles[TWIDSIZE];

int i;

/* Инициализация таблицы вращения */
twidfftf_fr32(twiddles,FFTSIZE);

/* Вычисление FFT от вещественного сигнала */
rfftf_fr32(sigdata, r_output, twiddles,1,FFTSIZE);

/* Добавление 2-й половины спектра */

for (i = 1; i < (FFTSIZE/2); i++)
{
   r_output[FFTSIZE-i] = conj_fr32(r_output[i]);
}

Функция вычисляет FFT над N вещественными выборками сигнала, используя алгоритм по основанию 4 (N-point radix-4 real input FFT).

#include < filter.h >

void rfftrad4_fr16(const fract16 input[],
                   complex_fract16 temp[],
                   complex_fract16 output[],
                   const complex_fract16 twiddle_table[],
                   int twiddle_stride,
                   int fft_size,
                   int block_exponent,
                   int scale_method);

Эта функция выполняет трансформацию входного сигнала домена времени в виде вещественных чисел в частотный домен с использованием БПФ по основанию 4. Достоинство функции rfftrad4_fr16 в том, что она использует факт равенства 0 мнимой части входного сигнала, что наполовину уменьшает количество умножений в алгоритме бабочки.

Размер входного буфера input, выходного буфера output и временного буфера temp равен fft_size, где fft_size представляет количество точек в FFT. Чтобы избежать потенциальных коллизий при обращениях к банку данных, буферы input и temp должны находиться в разных банках памяти. Это ускорит работу функции. Если входные данные могут быть перезаписаны, то оптимальное использование памяти может быть достигнуто также путем указания входного массива как выходного, при этом размер памяти входного массива должен быть как минимум 2*fft_size. Длина входной последовательности fft_size должна числом, равным степени двойки и не меньше 8.

Таблица вращения передается в аргументе twiddle_table, которая должна содержать как минимум 3*fft_size/4 множителей вращения. Для инициализации массива вращения может использоваться функция twidfftrad4_fr16. Если таблица вращения содержит больше множителей, чем нужно для конкретного вызова rfftrad4_fr16, то фактор шага таблицы вращения (параметр twiddle_stride) должен быть установлен соответствующим образом; иначе он должен быть установлен в 1.

Аргументы block_exponent и scale_method были добавлены с целью дальнейшего расширения функционала, в текущей реализации они игнорируются. Чтобы избежать переполнения, функция выполняет статическое масштабирование предварительным делением входа на fft_size.

Примечание: эта функция предоставлена для обратной совместимости с существующими приложениями. Новые приложения должны использовать функцию rfft_fr16 вместо функции rfftrad4_fr16.

[Алгоритм]

См. описание алгоритма функции cfftrad4.

Функция вычисляет магнитуду FFT.

#include < filter.h >

void fft_magnitude_fr16(const complex_fract16 input[],
                        fract16 output[],
                        int fft_size,
                        int block_exponent,
                        int mode);

void fft_magnitude_fr32(const complex_fract32 input[],
                        fract32 output[],
                        int fft_size,
                        int block_exponent,
                        int mode);

Функции магнитуды FFT fft_magnitude_fr16 и fft_magnitude_fr32 вычисляют нормализованный по мощности спектр (normalized power spectrum) от выходного сигнала, сгенерированного функцией FFT. аргумент fft_size задает размер FFT и должен быть числом, равным степени 2. Аргумент mode используется для указания типа функции FFT, которая использовалась для генерации входного массива. Функция fft_magnitude_fr16 вычисляет магнитуду FFT, которая представлена входным массивом с элементами типа fract16, в то время как fft_magnitude_fr32 вычисляет магнитуду FFT, которая представлена входным массивом с элементами типа fract32.

Если входной массив был сгенерирован из домена времени комплексного входного сигнала, параметр mode должен быть установлен в 0. Иначе аргумент mode должен быть установлен в 1, чтобы обозначить генерацию входного массива из домена времени вещественного входного сигнала. Т. е. mode=0 для входного массива, полученного вызовом библиотечных функций cfft_fr16, cfftf_fr16, cfftrad4_fr16, cfft_fr32, cfftf_fr32, и mode=1 для для входного массива, полученного вызовом функций rfft_fr16, rfftrad4_fr16, rfft_fr32, rfftf_fr32.

Аргумент block_exponent используется для управления нормализацией мощности спектра. Он обычно устанавливается в block_exponent, который возвращается функциями cfft_fr16 или cfft_fr32, rfft_fr16 или rfft_fr32. Если же входной массив был генерирован одной из функций cfftf_fr16 или cfftf_fr32, cfftrad4_fr16, rfftrad4_fr16 или rfftf_fr32, то аргумент block_exponent должен быть установлен в -1, показывая тем самым, что входной массив был генерирован с использованием статического масштабирования.

Если входной массив был сгенерирован какими-то другими средствами, то значение аргумента block_exponent будет зависеть от того, как было вычислено FFT. Если функция, использовавшаяся для вычисления FFT, не делала масштабирования промежуточных результатов на любом этапе вычисления, то установите block_exponent в 0; если функция FFT масштабировала результаты на каждом этапе вычислений, то установите block_exponent в -1; иначе установите block_exponent в количество стадий вычислений, которые делали масштабирование промежуточных результатов (это значение будет в диапазоне от 0 до log2(fft_size)).

Примечание: функции, которые вычисляют FFT с использованием арифметики с фиксированной точкой, будут обычно масштабировать набор промежуточных результатов для того, чтобы избежать для арифметики любых результатов с насыщением. См. описание функций cfft_fr16, rfft_fr16 или cfft_fr32, rfft_fr32 для получения дополнительной информации по различным методам масштабирования.

Функции fft_magnitude_fr16 и fft_magnitude_fr32 записывают мощность спектра в выходной массив output. Если mode установлено в 0, то длина мощности спектра будет fft_size. Если mode установлен в 1, то длина мощности спектра будет ((fft_size/2)+1).

[Условия возникновения ошибки]

Функции магнитуды FFT сделают немедленный выход без модификации вектора output, если верно любое из следующих условий:

• fft_size меньше 2.
• Аргумент mode установлен в значение, отличающееся от 0 или 1.
• block_exponent содержит значение меньше -1.
• block_exponent больше 0, и верно следующее: fft_size >= (1 << block_exponent).

[Алгоритм]

Для mode 0 (входной массив сгенерирован комплексными функциями cfft):

                     sqrt(input[i].re² + input[i].im²)
fft_magnitude[i] = --------------------------------------
                                 fft_size

где i = [0 ... fft_size).

Для mode 1 (входной массив сгенерирован вещественными функциями rfft):

                     2 * sqrt(input[i].re²)
fft_magnitude[i] = ---------------------------
                            fft_size

где i = [0 ... fft_size/2].

[Пример]

#include < filter.h >

#define N_FFT 1024

#pragma align 4096 

complex_fract16 cplx_signal[N_FFT];
fract16 real_signal[N_FFT];
complex_fract16 fft_output[N_FFT];
complex_fract16 twiddle_table[N_FFT];
fract16 real_magnitude[(N_FFT/2)+1];
fract16 cplx_magnitude[N_FFT];

int block_exponent; 

twidfftrad2_fr16(twiddle_table, N_FFT);
rfft_fr16(real_signal, fft_output, twiddle_table, 1, N_FFT, &block_exponent, 2);
fft_magnitude_fr16 (fft_output, real_magnitude, N_FFT, block_exponent, 1);
twidfftf_fr16 (twiddle_table, N_FFT);
cfftf_fr16(cplx_signal, fft_output, twiddle_table, 1, N_FFT);
fft_magnitude_fr16(fft_output, cplx_magnitude, N_FFT, -1, 0);

См. также: cfft, cfftf, rfft, rfftf.

Этот пример из проектов, которые поставляются вместе с VisualDSP++ (btc_fft_demo.c из проекта Background_Telemetry для платы разработчика ADSP-BF537 EZ-Kit Lite). Для диагностики и вывода результата используются так называемые каналы телеметрии [2] (Background Telemetry Channel, сокращенно BTC), поддержка которых встроена в VisualDSP++.

//--------------------------------------------------------------------------//

// Name: FFT BTC demo for the ADSP-BF537 EZ-KIT                             //

//--------------------------------------------------------------------------//

// Project Name: FFT Demo                                                   //

// Software: VisualDSP++ 4.0                                                //

// Hardware: ADSP-BF537 EZ-KIT Board                                        //

//--------------------------------------------------------------------------//

#include < sys/exception.h >

#include < complex.h >

#include < filter.h >

#include < math.h >

#include < math_const.h >

#include < btc.h >

#define NUMPOINTS 512

complex_fract16 out[NUMPOINTS];
complex_fract16 w[NUMPOINTS];       // Последовательность коэффициентов
                                    // косинусов и синусов (таблица вращения)

fract16 input_arr[NUMPOINTS];
fract16 mag[NUMPOINTS];

#define MINFREQ 100.0

#define MAXFREQ 20000.0

/////////////////////////

// Определения для BTC

/////////////////////////
fract16 BTC_CHAN0[NUMPOINTS+8];
fract16 BTC_CHAN1[NUMPOINTS+8];

int BTC_CHAN2 = 0x10;

BTC_MAP_BEGIN

//             имя канала, начальный адрес, длина
BTC_MAP_ENTRY("FFT_INPUT", (long)&BTC_CHAN0, sizeof(BTC_CHAN0))
BTC_MAP_ENTRY("FFT_OUTPUT", (long)&BTC_CHAN1, sizeof(BTC_CHAN1))
BTC_MAP_ENTRY("FREQ STEP SIZE", (long)&BTC_CHAN2, sizeof(BTC_CHAN2))
BTC_MAP_END

void initTimer(void);

void create_samples(float f);

EX_INTERRUPT_HANDLER(timerISR);        // обработчик прерывания таймера

void main()
{
   int i,j;
   int wst = 1;
   int n = NUMPOINTS;
   int block_exponent;
   int scale_method = 1;
   float freq;

   btc_init();
   // установка ISR
   register_handler(ik_timer, timerISR);
   // инициализация множителей вращения
   twidfftrad2_fr16(w, NUMPOINTS);
   //инкремент частоты
   freq = MINFREQ;
   create_samples(freq);
   // инициализация таймера и программируемых флагов
   initTimer();
   while (1)
   {
      // генерация входных выборок
      create_samples(freq);
      // БПФ массива выборок input_arr, спектр (в виде комплексных чисел)
      // будет помещен в массив out, используется таблица вращения w,
      // шаг по таблице wst=1, количество выборок n=512:
      rfft_fr16(input_arr, out, w, wst, n, &block_exponent, scale_method);
      for (i=0; i < NUMPOINTS; i++)
      {
         mag[i] = cabs_fr16(out[i]);
      }
      // запись в канал BTC:
      btc_write_array(0, (unsigned int*)input_arr, sizeof(input_arr));
      btc_write_array(1, (unsigned int*)mag, sizeof(mag));
      freq += BTC_CHAN2;
      if (freq > MAXFREQ)
         freq = MINFREQ;
   }
}

void create_samples(float f)
{
   int i,j;
   float fs = 48000.0;
   float step;
   //Генерация выборок для входных параметров БПФ:
   for (i=0; i < NUMPOINTS; i++)
   {
      step = (float)i/fs;
      input_arr[i] = (fract16) (sin(2*PI*f*step)*32760.0);
   }
}

void initTimer()
{
   unsigned int *mmrPtr;
   mmrPtr = (unsigned int*)0xffe03000;    // регистр управления таймера
   *mmrPtr = 5;
   mmrPtr = (unsigned int*)0xffe03004;    // регистр периода таймера
   *mmrPtr = 0x00001000;
   mmrPtr = (unsigned int*)0xffe03008;    // регистр масштаба таймера
   *mmrPtr = 0x00000000;
   mmrPtr = (unsigned int*)0xffe0300c;    // регистр счетчика таймера
   *mmrPtr = 0x00001000;
   mmrPtr = (unsigned int*)0xffe03000;    // регистр управления таймера
   *mmrPtr = 7;                           // разрешение работы таймера
}

EX_INTERRUPT_HANDLER(timerISR)
{
   btc_poll();
}

Этот пример выводит спектр сигнала до до 8 кГц, который оцифровывается с помощью драйвера SPORT (работа в режиме DMA) с частотой 128 кГц.

#include < filter.h >

//Константы для переменной mode, которые нужны для переключения режимов

// инициализации и запуска БПФ:

#define MODE_TEST_FFT_INIT    1

#define MODE_TEST_FFT_RUNNING 2

//Количество выборок сигнала в одном входном буфере:

#define ONEBUFSAMPLES 256

//Максимальное количество используемых входных буферов,

// от него зависят размеры массивов:

#define MAXNUMBUFFERS 1024

//Маска для указателей кольцевого буфера:

#define NUMBUFFMASK (MAXNUMBUFFERS-1)

//Функции БПФ будут использовать статическое масштабирование

// для промежуточных вычислений:

#define STATICSCALE 1

#define MINFREQ 20.0

#define MAXFREQ 8000.0

//Флажок, определяющий разрядность рассчитываемого БПФ:

// 16-битный БПФ или 32-битный БПФ:

bool fft16 = true;

//Массив для выходных данных, которые генерирует БПФ (спектр сигнала):

static section ("sdram0") complex_fract32 out32[MAXNUMBUFFERS*ONEBUFSAMPLES];

//Массив для таблицы вращения (используется функциями вычисления БПФ):

static section ("sdram0") complex_fract32 w32[MAXNUMBUFFERS*ONEBUFSAMPLES];

//Массив для входных данных БПФ (выборки сигнала):

static section ("sdram0") fract32 input_arr32[MAXNUMBUFFERS*ONEBUFSAMPLES];

//Указатели на аналогичные буферы, которые используются в 16-битном БПФ:
complex_fract16* out16;
complex_fract16* w16;
fract16* input_arr16;

//Шаг по частоте в спектре:

static float freqstep;

//Количество точек (выборок), по которым вычисляется БПФ:

static int numpoints;

//Количество буферов, которое используется для БПФ. По умолчанию

// 16 буферов, но можно задавать произвольное число, равное

// степени двойки, от 1 до MAXNUMBUFFERS.
u16 fftbuftotal = 16;

//Главный кольцевой буфер оцифровки, куда DMA льет данные

// из 18-битного АЦП блоками размером ONEBUFSAMPLES:

extern section ("sdram0") int samples[MAXNUMBUFFERS][ONEBUFSAMPLES];

extern u16 inSB, outSB;

//===========================================================================

// Функция, которая считает БПФ и отображает полученный линейный спектр

// на экране OLED. Эта функция вставляется либо в главный цикл main,

// либо в тело одного из потоков VDK (или другой RTOS):

//===========================================================================

void handlerFFTdemo (void)
{
   // Переменная, которая не используется. Сюда функция БПФ поместит
   // количество масштабирований, которое было произведено при вычислении
   // спектра:
   int block_exponent;
   
   if (MODE_TEST_FFT_INIT == mode)
   {

///////////////////////////////////////////////////////////////////

// Инициализация расчетов БПФ. В глобальную переменную fftbuftotal

// загружается количество используемых буферов, флажок fft16

// определяет, каким образом будет рассчитываться БПФ.
      numpoints = ONEBUFSAMPLES * fftbuftotal;
      // Вычисление шага спектра. 128 кГц частота оцифровки АЦП:
      freqstep = 128000.0/numpoints;
      // 16-битные буферы для экономии памяти используют область памяти,
      // выделенную под 32-битные буферы.
      out16       = (complex_fract16*)out32;
      w16         = (complex_fract16*)w32;
      input_arr16 = (fract16*)input_arr32;
      // Инициализация таблицы множителей вращения:
      if (fft16)
         twidfftrad2_fr16(w16, numpoints);
      else
         twidfftrad2_fr32(w32, numpoints);
      mode = MODE_TEST_FFT_RUNNING;
   }
   else
   {

///////////////////////////////////////////////////////////////////

// Расчет FFT по numpoints точкам.
      //Очистка экрана:
      u8g_SetColorIndex(0);
      ClearScreen();
      u8g_SetColorIndex(15);
LED(1);  //Миганием светодиод покажет начало и конец расчета
      // Смещение окна входных данных в на 1 буфер:
      u16 bufidxstart = outSB;
      outSB++;
      outSB &= NUMBUFFMASK;
      if (fft16)
      {
         // Подготовка 16-битного входного буфера для БПФ16:
         for (u16 bufcnt=0; bufcnt < fftbuftotal; bufcnt++)
         {
            for (int idx=0;idx < ONEBUFSAMPLES;idx++)
            {
               *(input_arr16+bufcnt*ONEBUFSAMPLES+idx) = 
                  (u32)samples[bufidxstart][idx]>>2;
            }
            bufidxstart++;
            bufidxstart &= NUMBUFFMASK;
         }
         // Вычисление БПФ16:
         rfft_fr16(input_arr16,     //входной массив
                   out16,           //выходной массив
                   w16,             //таблица вращения (cos и sin)
                   1,               //шаг по таблице вращения
                   numpoints,       //количество точек БПФ
                   &block_exponent, //вых. переменная масштабирования
                   STATICSCALE);    //используемое масштабирование для MAC
      }
      else
      {
         //Подготовка 32-битного входного буфера для БПФ32:
         for (u16 bufcnt=0; bufcnt < fftbuftotal; bufcnt++)
         {
            memcpy(&input_arr32[bufcnt*ONEBUFSAMPLES],
                   samples[bufidxstart],
                   ONEBUFSAMPLES*sizeof(int));
            bufidxstart++;
            bufidxstart &= NUMBUFFMASK;
         }
         // Вычисление БПФ32. Смысл параметров тот же:
         rfft_fr32(input_arr32,
                   out32,
                   w32,
                   1,
                   numpoints,
                   &block_exponent,
                   STATICSCALE);
      }
LED(0);  //Расчет спектра закончен

///////////////////////////////////////////////////////////////////

// Отображение спектра, который вычислило БПФ, на экране OLED.
      // Нулевую частоту спектра отбрасываем:
      int k = 1;
      // Цикл по 256 точкам оси X экрана:
      for (int x=0; x < 256; x++)
      {
         // Вычисление текущей частоты:
         float freq = k*freqstep;
         // Если частота больше заданной, то завершение:
         if (freq>MAXFREQ)
            break;
         // Рисование вертикальных палок в зависимости от амплитуды:
         // Функция cabs вычисляет абсолютное (вещественное) значение
         // амплитуды от комплексной величины спектра:
         int y;
         if(fft16)
            y = 63-cabs_fr16(out16[k])/32;
         else
            y = 63-cabs_fr32(out32[k])/128;
         if (y < 0)
            y=0;
         else if (y>63)
            y = 63;
         k++;
         // Вывод палки спектра на экран:
         u8g_DrawLine(x, 63, x, y);
      }
   }
}

Ниже показаны скриншотики спектра на экране OLED WEX025664.

На экране индикатора спектр от БПФ16 и БПФ32 выглядит совершенно одинаково, что свидетельствует о том, что расчеты выполняются правильно, и нет грубых ошибок в системе сбора и анализа данных. Однако время расчета БПФ16 и БПФ32 конечно отличается, хотя не нестолько, как я ожидал. Например, в режиме компиляции проекта Debug 16-битный БПФ по 16 буферам (4096 выборок сигнала) считается за 34 мс, в то время как 32-битный БПФ по тем же точкам считается 52 мс. В Release тот же БПФ вычисляется за 27.5 мс и 51 мс соответственно.