CMSIS DSP: функции БПФ для обработки действительных чисел Печать
Добавил(а) microsin   

В библиотеку CMSIS DSP включены специализированные алгоритмы для обработки FFT последовательности действительных чисел (Real FFT, RFFT). Преобразование FFT определено на комплексных числах, однако многие приложения в качестве входных данных имеют выборки действительных значений. Алгоритмы RFFT получают преимущество свойств симметрии FFT, и имеют преимущество над комплексными алгоритмами той же самой длины.

Быстрые алгоритмы RFFT полагаются на обработку смешанного основания (mixed-radix) CFFT, которая сохраняет процессорное время. На картинке ниже показаны шаги по вычислению прямого RFFT длины N.

[Real Fast Fourier Transform (RFFT)]

Действительная последовательность изначально обрабатывается как если бы она была комплексной, для выполнения CFFT. Далее на этапе обработки данные преобразуются для получения половины частотного спектра в комплексном формате. Кроме первого комплексного числа, которое содержит 2 действительных числа X[0] и X[N/2], все данные комплексные. Другими словами, первая комплексная выборка содержит 2 упакованных действительных значения.

RFFT

Рис. 1. Прямое преобразование БПФ для действительных чисел (Real Fast Fourier Transform, RFFT).

Вход инверсного RFFT должен иметь тот же формат, что и выходные данные прямого RFFT. На первом этапе обработки выполняется предварительная обработка данных для последующего выполнения обратного CFFT.

RIFFT

Рис. 2. Обратное преобразование БПФ для действительных чисел (Real Inverse Fast Fourier Transform, RIFFT).

Плавающая точка (float32). Основные функции arm_rfft_fast_f32() и arm_rfft_fast_init_f32(). Функции arm_rfft_f32() и arm_rfft_init_f32() устарели, но все еще документируются.

БПФ действительной N-точечной последовательности имеет четную симметрию в частотной области. Вторая половина данных равна сопряжению первой половины, перевернутой по частоте. При рассмотрении данных мы видим, что можно однозначно представлять FFT, используя только N/2 комплексных чисел. Они упакованы в выходном массиве, где чередуются действительные (real) и мнимые (imaginary) компоненты:

X = { real[0], imag[0], real[1], imag[1], real[2], imag[2] ... real[(N/2)-1], imag[(N/2)-1 }

Бывает, что первое комплексное значение (real[0], imag[0]) в действительности все real. real[0] представляет постоянное смещение, и imag[0] должен быть 0. (real[1], imag[1]) это фундаментальная частота, (real[2], imag[2]) первая гармоника и так далее.

Функции RFFT упаковывают данные домена частоты таким же способом. Прямое преобразование (из домена времени в домен частоты) выводит данные в этой форме, и обратное преобразование ожидает входные данные в этой же форме. Функция всегда выполняет необходимое реверсирование бит, чтобы входные и выходные данные всегда были в нормальном порядке. Функции поддерживают длины массивов [32, 64, 128, ..., 4096] выборок.

Фиксированная точка (Q15 и Q31). Алгоритмы RFFT определены аналогичным способом, и внутри себя они используют комплексное преобразование N/2. Комплексное преобразование внутри себя применяет масштабирование, чтобы избежать переполнений фиксированной точки. Общее масштабирование эквивалентно 1/(fftLen/2).

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

Существуют также соответствующие функции инициализации для каждого типа данных. Функция инициализации выполняет следующие действия:

• Устанавливает значения внутренних полей структуры.
• Инициализирует указатели на таблицы коэффициентов вращения и реверса бит.
• Инициализирует внутреннюю структуру данных комплексного FFT.

Использование функции инициализации не является обязательным. Однако, если она используется, то экземпляр структуры не может быть размещен в секцию данных const. Чтобы разместить экземпляр структуры в секцию данных const, то он должен быть инициализирован вручную следующим образом:

arm_rfft_instance_q31 S =
{
   fftLenReal,
   fftLenBy2,
   ifftFlagR,
   bitReverseFlagR,
   twidCoefRModifier,
   pTwiddleAReal,
   pTwiddleBReal,
   pCfft
};
 
arm_rfft_instance_q15 S =
{
   fftLenReal,
   fftLenBy2,
   ifftFlagR,
   bitReverseFlagR,
   twidCoefRModifier,
   pTwiddleAReal,
   pTwiddleBReal,
   pCfft
};

Здесь fftLenReal это длина преобразования последовательности действительных чисел, fftLenBy2 это длина внутреннего комплексного преобразования. ifftFlagR выбирает прямое (0) или инверсное (1) преобразование. bitReverseFlagR выбирает выход с реверсом бит (0) или нормальный порядок бит (1). twidCoefRModifier это модификатор шага для таблицы коэффициентов вращения. Значение основано на длине FFT. Указатель pTwiddleARealpoints указывает на массив A коэффициентов вращения; указатель pTwiddleBRealpoints указывает на массив B коэффициентов вращения; pCfft указывает на экземпляр структуры CFFT (Complex FFT). Структура CFFT также должна быть инициализирована. См. документацию arm_cfft_radix4_f32() для получения подробностей по статической инициализации экземпляра структуры Complex FFT.

[Описание функций RFFT CMSIS]

В библиотеку Real FFT CMSIS входят следующие функции:

// Обработка прямого или обратного преобразования Фурье
// (RFFT/RIFFT). Устаревший вариант.
void arm_rfft_f32 (const arm_rfft_instance_f32 *S,
                   float32_t *pSrc, float32_t *pDst);
 
// Прямое преобразование Фурье (RFFT) на числах с плавающей точкой.
void arm_rfft_fast_f32 (arm_rfft_fast_instance_f32 *S,
                        float32_t *p,
                        float32_t *pOut,
                        uint8_t ifftFlag);
 
// Функция инициализации для RFFT чисел с плавающей точкой.
arm_status arm_rfft_fast_init_f32 (arm_rfft_fast_instance_f32 *S,
                                   uint16_t fftLen);
 
// То же самое, но устаревший вариант.
arm_status arm_rfft_init_f32 (arm_rfft_instance_f32 *S,
                              arm_cfft_radix4_instance_f32 *S_CFFT,
                              uint32_t fftLenReal,
                              uint32_t ifftFlagR,
                              uint32_t bitReverseFlag);
 
// Функция инициализации для RFFT/RIFFT 16-битных
// чисел с фиксированной точкой (Q15).
arm_status arm_rfft_init_q15 (arm_rfft_instance_q15 *S,
                              uint32_t fftLenReal,
                              uint32_t ifftFlagR,
                              uint32_t bitReverseFlag);
 
// Функция инициализации для RFFT/RIFFT 32-битных
// чисел с фиксированной точкой (Q31).
arm_status arm_rfft_init_q31 (arm_rfft_instance_q31 *S,
                              uint32_t fftLenReal,
                              uint32_t ifftFlagR,
                              uint32_t bitReverseFlag);
 
// Функция прямого или обратного преобразования Фурье (RFFT/RIFFT)
// 16-битных чисел с фиксированной точкой (Q15).
void arm_rfft_q15 (const arm_rfft_instance_q15 *S,
                   q15_t *pSrc,
                   q15_t *pDst);
 
// Функция прямого или обратного преобразования Фурье (RFFT/RIFFT)
// 32-битных чисел с фиксированной точкой (Q31).
void arm_rfft_q31 (const arm_rfft_instance_q31 *S,
                   q31_t *pSrc, q31_t *pDst);

void arm_rfft_f32 (const arm_rfft_instance_f32 *S,
                   float32_t *pSrc, float32_t *pDst);

Прямое преобразование Фурье реальных чисел в формате с плавающей запятой. Устаревшая реализация функции, в будущих релизах библиотеки она будет удалена. Используйте вместо неё arm_rfft_fast_f32.

Параметры:

[in] *S указывает на (предварительно инициализированный) экземпляр структуры параметров RFFT/RIFFT.
[in] *pSrc указатель на входной буфер для обработки.
[out] *pDst указатель на выходной буфер.

Возвращаемое значение отсутствует.

void arm_rfft_fast_f32 (arm_rfft_fast_instance_f32 *S,
                        float32_t *p,
                        float32_t *pOut,
                        uint8_t ifftFlag);

Выполняет прямое или обратное преобразование Фурье для реальных чисел с плавающей запятой.

Параметры:

[in] *S указывает на (предварительно инициализированный) экземпляр структуры параметров RFFT/RIFFT типа arm_rfft_fast_instance_f32.
[in] *p указывает на входной буфер.
[in] *pOut указывает на выходной буфер.
[in] ifftFlag указывает, какое преобразование делать: 0 прямое (RFFT) или 1 обратное (RIFFT).

Примечание: прямое преобразование - из домена времени (выборки сигнала) в домен частот (спектр сигнала). Соответственно обратное обратное преобразование - из домена частот в домен сигнала.

arm_status arm_rfft_fast_init_f32 (arm_rfft_fast_instance_f32 *S,
                                   uint16_t fftLen);

Функция инициализирует структуру параметров для преобразования Фурье чисел с плавающей точкой.

Параметры:

[in, out] *S указывает на не иницилизированный экземпляр структуры arm_rfft_fast_instance_f32.
[in] fftLen длина обрабатываемого буфера (Real Sequence).

Вернет значение ARM_MATH_SUCCESS, если инициализация была успешна, или ARM_MATH_ARGUMENT_ERROR, если в параметре fftLen указано не поддерживаемое значение.

Параметр fftLen задает длину обработки RFFT/CIFFT. Поддерживаются фиксированные длины FFT 32, 64, 128, 256, 512, 1024, 2048, 4096. Функция также инициализирует указатель на таблицу коэффициентов вращения (Twiddle factor table) и указатель на таблицу реверса бит (Bit reversal table).

arm_status arm_rfft_init_f32 (arm_rfft_instance_f32 *S,
                              arm_cfft_radix4_instance_f32 *S_CFFT,
                              uint32_t fftLenReal,
                              uint32_t ifftFlagR,
                              uint32_t bitReverseFlag);

Эту функцию использует функция arm_dct4_init_f32().

arm_status arm_rfft_init_q15 (arm_rfft_instance_q15 *S,
                              uint32_t fftLenReal,
                              uint32_t ifftFlagR,
                              uint32_t bitReverseFlag);

Функция инициализации прямого или обратного преобразования Фурье для 15-битных чисел с фиксированной точкой.

Параметры:

[in, out] *S указывает на не иницилизированный экземпляр структуры Q15 RFFT/RIFFT.
[in] fftLenReal длина FFT.
[in] ifftFlagR указывает, какое преобразование делать: 0 прямое (RFFT) или 1 обратное (RIFFT).
[in] bitReverseFlag флаг, который разрешает инверсию порядка следования бит на выходе (bitReverseFlag=1) или запрещает её (bitReverseFlag=0).

Вернет значение ARM_MATH_SUCCESS, если инициализация была успешна, или ARM_MATH_ARGUMENT_ERROR, если в параметре fftLen указано не поддерживаемое значение.

Параметр fftLenReal задает длину обработки RFFT/RIFFT. Поддерживаются фиксированные длины FFT 32, 64, 128, 256, 512, 1024, 2048, 4096, 8192. Функция также инициализирует таблицу коэффициентов вращения (Twiddle factor table) и таблицу реверса бит (Bit reversal table).

arm_status arm_rfft_init_q31 (arm_rfft_instance_q31 *S,
                              uint32_t fftLenReal,
                              uint32_t ifftFlagR,
                              uint32_t bitReverseFlag);

То же самое, что и arm_rfft_init_q15, только для 32-битного формата чисел с фиксированной запятой Q31.

void arm_rfft_q15 (const arm_rfft_instance_q15 *S,
                   q15_t *pSrc,
                   q15_t *pDst);

Функция выполняет прямое или обратное преобразование Фурье для 16-битных чисел с фиксированной запятой. Следует заметить, что по обработки она не лучше функции с плавающей точкой arm_rfft_fast_f32. Однако у arm_rfft_init_q15 есть преимущество в том, что поддерживаются буферы размера 8192.

Параметры:

[in] *S указывает на (предварительно инициализированный) экземпляр структуры Q15 RFFT/RIFFT.
[in] *pSrc указывает на входной буфер.
[out] *pDst указывает на выходной буфер.

Возвращаемое значение отсутствует.

[Входные и выходные форматы данных]

Внутри функции данные масштабируются в сторону уменьшения в 2 раза на каждом шаге, что делается с целью избежать переполнения (насыщения) чисел с фиксированной запятой в процессе CFFT/CIFFT. Поэтому представление данных на выходе отличается для разных размеров RFFT. Формат на входе и выходе различных размеров RFFT и количество бит при масштабировании RFFT и RIFFT показаны в таблицах ниже.

Входные и выходные форматы для Q15 RFFT:

Размер Входной формат Выходной формат Сколько бит теряется
32 1.15 5.11 4
64 6.10 5
128 7.9 6
256 8.8 7
512 9.7 8
1024 10.6 9
2048 11.5 10
4096 12.4 11
8192 13.3 12

Входные и выходные форматы для Q15 RIFFT:

Размер Входной формат Выходной формат Сколько бит теряется
32 1.15 5.11 0
64 6.10 0
128 7.9 0
256 8.8 0
512 9.7 0
1024 10.6 0
2048 11.5 0
4096 12.4 0
8192 13.3 0

[Пример использования arm_rfft_q15]

#include < math.h>
#include < string.h>
 
#define FFT_LENGTH 8192
 
short rawADCbuffer[FFT_LENGTH];        // Исходные сырые данные звука
static q15_t bufFFTi[FFT_LENGTH*2];    // Входной буфер FFT (сигнал)
static q15_t bufFFTo[FFT_LENGTH*2];    // Выходной буфер FFT (спектр)
static float maxMean = 0.;
static float maxValue = 0.;
 
...
// Подготовка буфера q15_t из данных АЦП в формате short (int16_t):
for (int i=0; i < FFT_LENGTH; i++)
{
   bufFFTi[i] = (q15_t)(rawADCbuffer[i]);
   bufFFTi[i+FFT_LENGTH] = 0;
}
 
// Подготовка структуры параметров для RFFT.
arm_rfft_instance_q15 rffti;
// 0:прямое преобразование (forward FT), 1: нормальный порядок следования бит
arm_rfft_init_q15(&rffti, FFT_LENGTH, 0, 1);
 
// Выполнение преобразования Фурье:
arm_rfft_q15(&rffti, bufFFTi, bufFFTo);
 
/* Обработка данных с помощью Complex Magnitude Module, чтобы
   вычислить магнитуду каждого бина. Для экономии памяти bufFFTo
   теперь служит входным буфером, где находятся данные спектра,
   а в bufFFTi будут находиться реальные данные магнитуд. */
arm_cmplx_mag_q15(bufFFTo, bufFFTi, FFT_LENGTH);
 
/* Вычисление maxValue и возврат значения соответствующего бина: */
q15_t maxValue_q15 = 0;uint32_t maxIndex = 0;
maxValue = 0.;
// Здесь сканируется не весь массив спектра, а только
// интересующий кусок от индекса частоты filterFreqIndexLow
// до индекса частоты filterFreqIndexHigh. Максимальное значение
// сохраняется в maxValue_q15, а его индекс частоты в maxIndex.
arm_max_q15(bufFFTi + filterFreqIndexLow,
            filterFreqIndexHigh - filterFreqIndexLow,
            &maxValue_q15,
            &maxIndex);
maxIndex += filterFreqIndexLow;
 
// Преобразование q15 во float:
arm_q15_to_float(&maxValue_q15, &maxValue, 1);
maxValue *= 10;
 
// Вычисление среднего значения с добавлением соседних бинов:
float maxValuePrev = 0.;
arm_q15_to_float(&bufFFTi[maxIndex-1], &maxValuePrev, 1);
maxValuePrev *= 10;float maxValueNext = 0.;
arm_q15_to_float(&bufFFTi[maxIndex+1], &maxValueNext, 1);
maxValueNext *= 10;if(maxValue > 0.05)
   maxMean += (maxValueNext - maxValuePrev) / maxValue;

void arm_rfft_q31 (const arm_rfft_instance_q31 *S,
                   q31_t *pSrc, q31_t *pDst);

Функция выполняет прямое или обратное преобразование Фурье для 32-битных чисел с фиксированной запятой. К сожалению, на STM32 скорость работы этой функции оставляет желать лучшего, и она не подойдет для обработки данных на лету, в реальном времени.

Параметры:

[in] *S указывает на (предварительно инициализированный) экземпляр структуры Q31 RFFT/RIFFT.
[in] *pSrc указывает на входной буфер.
[out] *pDst указывает на выходной буфер.

Возвращаемое значение отсутствует.

[Входные и выходные форматы данных]

Внутри функции данные масштабируются в сторону уменьшения в 2 раза на каждом шаге, что делается с целью избежать переполнения (насыщения) чисел с фиксированной запятой в процессе CFFT/CIFFT. Поэтому представление данных на выходе отличается для разных размеров RFFT. Формат на входе и выходе различных размеров RFFT и количество бит при масштабировании RFFT и RIFFT показаны в таблицах ниже.

Входные и выходные форматы для Q31 RFFT:

Размер Входной формат Выходной формат Сколько бит теряется
32 1.31 5.27 4
64 6.26 5
128 7.25 6
256 8.24 7
512 9.23 8
1024 10.22 9
2048 11.21 10
4096 12.20 11
8192 13.19 12

Входные и выходные форматы для Q31 RIFFT:

Размер Входной формат Выходной формат Сколько бит теряется
32 1.31 5.27 0
64 6.26 0
128 7.25 0
256 8.24 0
512 9.23 0
1024 10.22 0
2048 11.21 0
4096 12.20 0
8192 13.19 0

[Таблицы Real FFT]

Переменные таблиц:

static const float32_t realCoefA [8192];
static const float32_t realCoefB [8192];
const q15_t ALIGN4 realCoefAQ15 [8192];
const q15_t ALIGN4 realCoefBQ15 [8192];
const q31_t realCoefAQ31 [8192];
const q31_t realCoefBQ31 [8192];

Таблицы realCoefA, realCoefB используются устаревшей функцией arm_rfft_init_f32().

Алгоритм генерации realCoefA:

n = 4096;for (i = 0; i < n; i++)
{
   pATable[2 * i] = 0.5 * (1.0 - sin (2 * PI / (double) (2 * n) * (double) i));
   pATable[2 * i + 1] = 0.5 * (-1.0 * cos (2 * PI / (double) (2 * n) * (double) i));
}

Алгоритм генерации realCoefB:

n = 4096;for (i = 0; i < n; i++)
{
   pBTable[2 * i] = 0.5 * (1.0 + sin (2 * PI / (double) (2 * n) * (double) i));
   pBTable[2 * i + 1] = 0.5 * (1.0 * cos (2 * PI / (double) (2 * n) * (double) i));
}

Таблицы realCoefAQ15 и realCoefBQ15 используются функцией arm_rfft_init_q15(). Таблицы realCoefAQ31 и realCoefBQ31 используются функцией arm_rfft_init_q31(). Алгоритм генерации этих таблиц такой же, только коэффициенты преобразуются в целочисленный формат операциями round(pATable[i] * pow(2, 15)) и round(floatТаблица[i] * pow(2, 31)) для форматов Q15 и Q31 соответственно.

[Ссылки]

1. Real FFT Functions site:keil.com.
2. CMSIS DSP: комплексные функции БПФ.