luv090909/py32f0-template
1
0
Fork 0
mirror of https://github.com/IcedRooibos/py32f0-template.git synced 2025-10-31 01:42:04 -07:00
Code Issues Packages Projects Releases Wiki Activity
py32f0-template/Libraries/CMSIS/DSP/Source/TransformFunctions/arm_cfft_f32.c
IOsetting e19b6df4e3 fix: incorrect file attributes
2023-01-19 22:45:57 +08:00

1193 lines
35 KiB
C
Raw Blame History

/* ----------------------------------------------------------------------
* Project: CMSIS DSP Library
* Title: arm_cfft_f32.c
* Description: Combined Radix Decimation in Frequency CFFT Floating point processing function
*
* $Date: 23 April 2021
* $Revision: V1.9.0
*
* Target Processor: Cortex-M and Cortex-A cores
* -------------------------------------------------------------------- */
/*
* Copyright (C) 2010-2021 ARM Limited or its affiliates. All rights reserved.
*
* SPDX-License-Identifier: Apache-2.0
*
* Licensed under the Apache License, Version 2.0 (the License); you may
* not use this file except in compliance with the License.
* You may obtain a copy of the License at
*
* www.apache.org/licenses/LICENSE-2.0
*
* Unless required by applicable law or agreed to in writing, software
* distributed under the License is distributed on an AS IS BASIS, WITHOUT
* WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
* See the License for the specific language governing permissions and
* limitations under the License.
*/
#include "dsp/transform_functions.h"
#include "arm_common_tables.h"
#if defined(ARM_MATH_MVEF) && !defined(ARM_MATH_AUTOVECTORIZE)
#include "arm_helium_utils.h"
#include "arm_vec_fft.h"
#include "arm_mve_tables.h"
static float32_t arm_inverse_fft_length_f32(uint16_t fftLen)
{
float32_t retValue=1.0;
switch (fftLen)
{
case 4096U:
retValue = 0.000244140625;
break;
case 2048U:
retValue = 0.00048828125;
break;
case 1024U:
retValue = 0.0009765625f;
break;
case 512U:
retValue = 0.001953125;
break;
case 256U:
retValue = 0.00390625f;
break;
case 128U:
retValue = 0.0078125;
break;
case 64U:
retValue = 0.015625f;
break;
case 32U:
retValue = 0.03125;
break;
case 16U:
retValue = 0.0625f;
break;
default:
break;
}
return(retValue);
}
static void _arm_radix4_butterfly_f32_mve(const arm_cfft_instance_f32 * S,float32_t * pSrc, uint32_t fftLen)
{
f32x4_t vecTmp0, vecTmp1;
f32x4_t vecSum0, vecDiff0, vecSum1, vecDiff1;
f32x4_t vecA, vecB, vecC, vecD;
uint32_t blkCnt;
uint32_t n1, n2;
uint32_t stage = 0;
int32_t iter = 1;
static const int32_t strides[4] = {
(0 - 16) * (int32_t)sizeof(q31_t *),
(1 - 16) * (int32_t)sizeof(q31_t *),
(8 - 16) * (int32_t)sizeof(q31_t *),
(9 - 16) * (int32_t)sizeof(q31_t *)
};
n2 = fftLen;
n1 = n2;
n2 >>= 2u;
for (int k = fftLen / 4u; k > 1; k >>= 2)
{
float32_t const *p_rearranged_twiddle_tab_stride1 =
&S->rearranged_twiddle_stride1[
S->rearranged_twiddle_tab_stride1_arr[stage]];
float32_t const *p_rearranged_twiddle_tab_stride2 =
&S->rearranged_twiddle_stride2[
S->rearranged_twiddle_tab_stride2_arr[stage]];
float32_t const *p_rearranged_twiddle_tab_stride3 =
&S->rearranged_twiddle_stride3[
S->rearranged_twiddle_tab_stride3_arr[stage]];
float32_t * pBase = pSrc;
for (int i = 0; i < iter; i++)
{
float32_t *inA = pBase;
float32_t *inB = inA + n2 * CMPLX_DIM;
float32_t *inC = inB + n2 * CMPLX_DIM;
float32_t *inD = inC + n2 * CMPLX_DIM;
float32_t const *pW1 = p_rearranged_twiddle_tab_stride1;
float32_t const *pW2 = p_rearranged_twiddle_tab_stride2;
float32_t const *pW3 = p_rearranged_twiddle_tab_stride3;
f32x4_t vecW;
blkCnt = n2 / 2;
/*
* load 2 f32 complex pair
*/
vecA = vldrwq_f32(inA);
vecC = vldrwq_f32(inC);
while (blkCnt > 0U)
{
vecB = vldrwq_f32(inB);
vecD = vldrwq_f32(inD);
vecSum0 = vecA + vecC; /* vecSum0 = vaddq(vecA, vecC) */
vecDiff0 = vecA - vecC; /* vecSum0 = vsubq(vecA, vecC) */
vecSum1 = vecB + vecD;
vecDiff1 = vecB - vecD;
/*
* [ 1 1 1 1 ] * [ A B C D ]' .* 1
*/
vecTmp0 = vecSum0 + vecSum1;
vst1q(inA, vecTmp0);
inA += 4;
/*
* [ 1 -1 1 -1 ] * [ A B C D ]'
*/
vecTmp0 = vecSum0 - vecSum1;
/*
* [ 1 -1 1 -1 ] * [ A B C D ]'.* W2
*/
vecW = vld1q(pW2);
pW2 += 4;
vecTmp1 = MVE_CMPLX_MULT_FLT_Conj_AxB(vecW, vecTmp0);
vst1q(inB, vecTmp1);
inB += 4;
/*
* [ 1 -i -1 +i ] * [ A B C D ]'
*/
vecTmp0 = MVE_CMPLX_SUB_A_ixB(vecDiff0, vecDiff1);
/*
* [ 1 -i -1 +i ] * [ A B C D ]'.* W1
*/
vecW = vld1q(pW1);
pW1 +=4;
vecTmp1 = MVE_CMPLX_MULT_FLT_Conj_AxB(vecW, vecTmp0);
vst1q(inC, vecTmp1);
inC += 4;
/*
* [ 1 +i -1 -i ] * [ A B C D ]'
*/
vecTmp0 = MVE_CMPLX_ADD_A_ixB(vecDiff0, vecDiff1);
/*
* [ 1 +i -1 -i ] * [ A B C D ]'.* W3
*/
vecW = vld1q(pW3);
pW3 += 4;
vecTmp1 = MVE_CMPLX_MULT_FLT_Conj_AxB(vecW, vecTmp0);
vst1q(inD, vecTmp1);
inD += 4;
vecA = vldrwq_f32(inA);
vecC = vldrwq_f32(inC);
blkCnt--;
}
pBase += CMPLX_DIM * n1;
}
n1 = n2;
n2 >>= 2u;
iter = iter << 2;
stage++;
}
/*
* start of Last stage process
*/
uint32x4_t vecScGathAddr = vld1q_u32((uint32_t*)strides);
vecScGathAddr = vecScGathAddr + (uint32_t) pSrc;
/* load scheduling */
vecA = vldrwq_gather_base_wb_f32(&vecScGathAddr, 64);
vecC = vldrwq_gather_base_f32(vecScGathAddr, 16);
blkCnt = (fftLen >> 3);
while (blkCnt > 0U)
{
vecSum0 = vecA + vecC; /* vecSum0 = vaddq(vecA, vecC) */
vecDiff0 = vecA - vecC; /* vecSum0 = vsubq(vecA, vecC) */
vecB = vldrwq_gather_base_f32(vecScGathAddr, 8);
vecD = vldrwq_gather_base_f32(vecScGathAddr, 24);
vecSum1 = vecB + vecD;
vecDiff1 = vecB - vecD;
/* pre-load for next iteration */
vecA = vldrwq_gather_base_wb_f32(&vecScGathAddr, 64);
vecC = vldrwq_gather_base_f32(vecScGathAddr, 16);
vecTmp0 = vecSum0 + vecSum1;
vstrwq_scatter_base_f32(vecScGathAddr, -64, vecTmp0);
vecTmp0 = vecSum0 - vecSum1;
vstrwq_scatter_base_f32(vecScGathAddr, -64 + 8, vecTmp0);
vecTmp0 = MVE_CMPLX_SUB_A_ixB(vecDiff0, vecDiff1);
vstrwq_scatter_base_f32(vecScGathAddr, -64 + 16, vecTmp0);
vecTmp0 = MVE_CMPLX_ADD_A_ixB(vecDiff0, vecDiff1);
vstrwq_scatter_base_f32(vecScGathAddr, -64 + 24, vecTmp0);
blkCnt--;
}
/*
* End of last stage process
*/
}
static void arm_cfft_radix4by2_f32_mve(const arm_cfft_instance_f32 * S, float32_t *pSrc, uint32_t fftLen)
{
float32_t const *pCoefVec;
float32_t const *pCoef = S->pTwiddle;
float32_t *pIn0, *pIn1;
uint32_t n2;
uint32_t blkCnt;
f32x4_t vecIn0, vecIn1, vecSum, vecDiff;
f32x4_t vecCmplxTmp, vecTw;
n2 = fftLen >> 1;
pIn0 = pSrc;
pIn1 = pSrc + fftLen;
pCoefVec = pCoef;
blkCnt = n2 / 2;
while (blkCnt > 0U)
{
vecIn0 = *(f32x4_t *) pIn0;
vecIn1 = *(f32x4_t *) pIn1;
vecTw = vld1q(pCoefVec);
pCoefVec += 4;
vecSum = vecIn0 + vecIn1;
vecDiff = vecIn0 - vecIn1;
vecCmplxTmp = MVE_CMPLX_MULT_FLT_Conj_AxB(vecTw, vecDiff);
vst1q(pIn0, vecSum);
pIn0 += 4;
vst1q(pIn1, vecCmplxTmp);
pIn1 += 4;
blkCnt--;
}
_arm_radix4_butterfly_f32_mve(S, pSrc, n2);
_arm_radix4_butterfly_f32_mve(S, pSrc + fftLen, n2);
pIn0 = pSrc;
}
static void _arm_radix4_butterfly_inverse_f32_mve(const arm_cfft_instance_f32 * S,float32_t * pSrc, uint32_t fftLen, float32_t onebyfftLen)
{
f32x4_t vecTmp0, vecTmp1;
f32x4_t vecSum0, vecDiff0, vecSum1, vecDiff1;
f32x4_t vecA, vecB, vecC, vecD;
uint32_t blkCnt;
uint32_t n1, n2;
uint32_t stage = 0;
int32_t iter = 1;
static const int32_t strides[4] = {
(0 - 16) * (int32_t)sizeof(q31_t *),
(1 - 16) * (int32_t)sizeof(q31_t *),
(8 - 16) * (int32_t)sizeof(q31_t *),
(9 - 16) * (int32_t)sizeof(q31_t *)
};
n2 = fftLen;
n1 = n2;
n2 >>= 2u;
for (int k = fftLen / 4; k > 1; k >>= 2)
{
float32_t const *p_rearranged_twiddle_tab_stride1 =
&S->rearranged_twiddle_stride1[
S->rearranged_twiddle_tab_stride1_arr[stage]];
float32_t const *p_rearranged_twiddle_tab_stride2 =
&S->rearranged_twiddle_stride2[
S->rearranged_twiddle_tab_stride2_arr[stage]];
float32_t const *p_rearranged_twiddle_tab_stride3 =
&S->rearranged_twiddle_stride3[
S->rearranged_twiddle_tab_stride3_arr[stage]];
float32_t * pBase = pSrc;
for (int i = 0; i < iter; i++)
{
float32_t *inA = pBase;
float32_t *inB = inA + n2 * CMPLX_DIM;
float32_t *inC = inB + n2 * CMPLX_DIM;
float32_t *inD = inC + n2 * CMPLX_DIM;
float32_t const *pW1 = p_rearranged_twiddle_tab_stride1;
float32_t const *pW2 = p_rearranged_twiddle_tab_stride2;
float32_t const *pW3 = p_rearranged_twiddle_tab_stride3;
f32x4_t vecW;
blkCnt = n2 / 2;
/*
* load 2 f32 complex pair
*/
vecA = vldrwq_f32(inA);
vecC = vldrwq_f32(inC);
while (blkCnt > 0U)
{
vecB = vldrwq_f32(inB);
vecD = vldrwq_f32(inD);
vecSum0 = vecA + vecC; /* vecSum0 = vaddq(vecA, vecC) */
vecDiff0 = vecA - vecC; /* vecSum0 = vsubq(vecA, vecC) */
vecSum1 = vecB + vecD;
vecDiff1 = vecB - vecD;
/*
* [ 1 1 1 1 ] * [ A B C D ]' .* 1
*/
vecTmp0 = vecSum0 + vecSum1;
vst1q(inA, vecTmp0);
inA += 4;
/*
* [ 1 -1 1 -1 ] * [ A B C D ]'
*/
vecTmp0 = vecSum0 - vecSum1;
/*
* [ 1 -1 1 -1 ] * [ A B C D ]'.* W1
*/
vecW = vld1q(pW2);
pW2 += 4;
vecTmp1 = MVE_CMPLX_MULT_FLT_AxB(vecW, vecTmp0);
vst1q(inB, vecTmp1);
inB += 4;
/*
* [ 1 -i -1 +i ] * [ A B C D ]'
*/
vecTmp0 = MVE_CMPLX_ADD_A_ixB(vecDiff0, vecDiff1);
/*
* [ 1 -i -1 +i ] * [ A B C D ]'.* W2
*/
vecW = vld1q(pW1);
pW1 += 4;
vecTmp1 = MVE_CMPLX_MULT_FLT_AxB(vecW, vecTmp0);
vst1q(inC, vecTmp1);
inC += 4;
/*
* [ 1 +i -1 -i ] * [ A B C D ]'
*/
vecTmp0 = MVE_CMPLX_SUB_A_ixB(vecDiff0, vecDiff1);
/*
* [ 1 +i -1 -i ] * [ A B C D ]'.* W3
*/
vecW = vld1q(pW3);
pW3 += 4;
vecTmp1 = MVE_CMPLX_MULT_FLT_AxB(vecW, vecTmp0);
vst1q(inD, vecTmp1);
inD += 4;
vecA = vldrwq_f32(inA);
vecC = vldrwq_f32(inC);
blkCnt--;
}
pBase += CMPLX_DIM * n1;
}
n1 = n2;
n2 >>= 2u;
iter = iter << 2;
stage++;
}
/*
* start of Last stage process
*/
uint32x4_t vecScGathAddr = vld1q_u32 ((uint32_t*)strides);
vecScGathAddr = vecScGathAddr + (uint32_t) pSrc;
/*
* load scheduling
*/
vecA = vldrwq_gather_base_wb_f32(&vecScGathAddr, 64);
vecC = vldrwq_gather_base_f32(vecScGathAddr, 16);
blkCnt = (fftLen >> 3);
while (blkCnt > 0U)
{
vecSum0 = vecA + vecC; /* vecSum0 = vaddq(vecA, vecC) */
vecDiff0 = vecA - vecC; /* vecSum0 = vsubq(vecA, vecC) */
vecB = vldrwq_gather_base_f32(vecScGathAddr, 8);
vecD = vldrwq_gather_base_f32(vecScGathAddr, 24);
vecSum1 = vecB + vecD;
vecDiff1 = vecB - vecD;
vecA = vldrwq_gather_base_wb_f32(&vecScGathAddr, 64);
vecC = vldrwq_gather_base_f32(vecScGathAddr, 16);
vecTmp0 = vecSum0 + vecSum1;
vecTmp0 = vecTmp0 * onebyfftLen;
vstrwq_scatter_base_f32(vecScGathAddr, -64, vecTmp0);
vecTmp0 = vecSum0 - vecSum1;
vecTmp0 = vecTmp0 * onebyfftLen;
vstrwq_scatter_base_f32(vecScGathAddr, -64 + 8, vecTmp0);
vecTmp0 = MVE_CMPLX_ADD_A_ixB(vecDiff0, vecDiff1);
vecTmp0 = vecTmp0 * onebyfftLen;
vstrwq_scatter_base_f32(vecScGathAddr, -64 + 16, vecTmp0);
vecTmp0 = MVE_CMPLX_SUB_A_ixB(vecDiff0, vecDiff1);
vecTmp0 = vecTmp0 * onebyfftLen;
vstrwq_scatter_base_f32(vecScGathAddr, -64 + 24, vecTmp0);
blkCnt--;
}
/*
* End of last stage process
*/
}
static void arm_cfft_radix4by2_inverse_f32_mve(const arm_cfft_instance_f32 * S,float32_t *pSrc, uint32_t fftLen)
{
float32_t const *pCoefVec;
float32_t const *pCoef = S->pTwiddle;
float32_t *pIn0, *pIn1;
uint32_t n2;
float32_t onebyfftLen = arm_inverse_fft_length_f32(fftLen);
uint32_t blkCnt;
f32x4_t vecIn0, vecIn1, vecSum, vecDiff;
f32x4_t vecCmplxTmp, vecTw;
n2 = fftLen >> 1;
pIn0 = pSrc;
pIn1 = pSrc + fftLen;
pCoefVec = pCoef;
blkCnt = n2 / 2;
while (blkCnt > 0U)
{
vecIn0 = *(f32x4_t *) pIn0;
vecIn1 = *(f32x4_t *) pIn1;
vecTw = vld1q(pCoefVec);
pCoefVec += 4;
vecSum = vecIn0 + vecIn1;
vecDiff = vecIn0 - vecIn1;
vecCmplxTmp = MVE_CMPLX_MULT_FLT_AxB(vecTw, vecDiff);
vst1q(pIn0, vecSum);
pIn0 += 4;
vst1q(pIn1, vecCmplxTmp);
pIn1 += 4;
blkCnt--;
}
_arm_radix4_butterfly_inverse_f32_mve(S, pSrc, n2, onebyfftLen);
_arm_radix4_butterfly_inverse_f32_mve(S, pSrc + fftLen, n2, onebyfftLen);
}
/**
@addtogroup ComplexFFT
@{
*/
/**
@brief Processing function for the floating-point complex FFT.
@param[in] S points to an instance of the floating-point CFFT structure
@param[in,out] p1 points to the complex data buffer of size <code>2*fftLen</code>. Processing occurs in-place
@param[in] ifftFlag flag that selects transform direction
- value = 0: forward transform
- value = 1: inverse transform
@param[in] bitReverseFlag flag that enables / disables bit reversal of output
- value = 0: disables bit reversal of output
- value = 1: enables bit reversal of output
@return none
*/
void arm_cfft_f32(
const arm_cfft_instance_f32 * S,
float32_t * pSrc,
uint8_t ifftFlag,
uint8_t bitReverseFlag)
{
uint32_t fftLen = S->fftLen;
if (ifftFlag == 1U) {
switch (fftLen) {
case 16:
case 64:
case 256:
case 1024:
case 4096:
_arm_radix4_butterfly_inverse_f32_mve(S, pSrc, fftLen, arm_inverse_fft_length_f32(S->fftLen));
break;
case 32:
case 128:
case 512:
case 2048:
arm_cfft_radix4by2_inverse_f32_mve(S, pSrc, fftLen);
break;
}
} else {
switch (fftLen) {
case 16:
case 64:
case 256:
case 1024:
case 4096:
_arm_radix4_butterfly_f32_mve(S, pSrc, fftLen);
break;
case 32:
case 128:
case 512:
case 2048:
arm_cfft_radix4by2_f32_mve(S, pSrc, fftLen);
break;
}
}
if (bitReverseFlag)
{
arm_bitreversal_32_inpl_mve((uint32_t*)pSrc, S->bitRevLength, S->pBitRevTable);
}
}
#else
extern void arm_radix8_butterfly_f32(
float32_t * pSrc,
uint16_t fftLen,
const float32_t * pCoef,
uint16_t twidCoefModifier);
extern void arm_bitreversal_32(
uint32_t * pSrc,
const uint16_t bitRevLen,
const uint16_t * pBitRevTable);
/**
@ingroup groupTransforms
*/
/**
@defgroup ComplexFFT Complex FFT Functions
@par
The Fast Fourier Transform (FFT) is an efficient algorithm for computing the
Discrete Fourier Transform (DFT). The FFT can be orders of magnitude faster
than the DFT, especially for long lengths.
The algorithms described in this section
operate on complex data. A separate set of functions is devoted to handling
of real sequences.
@par
There are separate algorithms for handling floating-point, Q15, and Q31 data
types. The algorithms available for each data type are described next.
@par
The FFT functions operate in-place. That is, the array holding the input data
will also be used to hold the corresponding result. The input data is complex
and contains <code>2*fftLen</code> interleaved values as shown below.
<pre>{real[0], imag[0], real[1], imag[1], ...} </pre>
The FFT result will be contained in the same array and the frequency domain
values will have the same interleaving.
@par Floating-point
The floating-point complex FFT uses a mixed-radix algorithm. Multiple radix-8
stages are performed along with a single radix-2 or radix-4 stage, as needed.
The algorithm supports lengths of [16, 32, 64, ..., 4096] and each length uses
a different twiddle factor table.
@par
The function uses the standard FFT definition and output values may grow by a
factor of <code>fftLen</code> when computing the forward transform. The
inverse transform includes a scale of <code>1/fftLen</code> as part of the
calculation and this matches the textbook definition of the inverse FFT.
@par
For the MVE version, the new arm_cfft_init_f32 initialization function is
<b>mandatory</b>. <b>Compilation flags are available to include only the required tables for the
needed FFTs.</b> Other FFT versions can continue to be initialized as
explained below.
@par
For not MVE versions, pre-initialized data structures containing twiddle factors
and bit reversal tables are provided and defined in <code>arm_const_structs.h</code>. Include
this header in your function and then pass one of the constant structures as
an argument to arm_cfft_f32. For example:
@par
<code>arm_cfft_f32(arm_cfft_sR_f32_len64, pSrc, 1, 1)</code>
@par
computes a 64-point inverse complex FFT including bit reversal.
The data structures are treated as constant data and not modified during the
calculation. The same data structure can be reused for multiple transforms
including mixing forward and inverse transforms.
@par
Earlier releases of the library provided separate radix-2 and radix-4
algorithms that operated on floating-point data. These functions are still
provided but are deprecated. The older functions are slower and less general
than the new functions.
@par
An example of initialization of the constants for the arm_cfft_f32 function follows:
@code
const static arm_cfft_instance_f32 *S;
...
switch (length) {
case 16:
S = &arm_cfft_sR_f32_len16;
break;
case 32:
S = &arm_cfft_sR_f32_len32;
break;
case 64:
S = &arm_cfft_sR_f32_len64;
break;
case 128:
S = &arm_cfft_sR_f32_len128;
break;
case 256:
S = &arm_cfft_sR_f32_len256;
break;
case 512:
S = &arm_cfft_sR_f32_len512;
break;
case 1024:
S = &arm_cfft_sR_f32_len1024;
break;
case 2048:
S = &arm_cfft_sR_f32_len2048;
break;
case 4096:
S = &arm_cfft_sR_f32_len4096;
break;
}
@endcode
@par
The new arm_cfft_init_f32 can also be used.
@par Q15 and Q31
The floating-point complex FFT uses a mixed-radix algorithm. Multiple radix-4
stages are performed along with a single radix-2 stage, as needed.
The algorithm supports lengths of [16, 32, 64, ..., 4096] and each length uses
a different twiddle factor table.
@par
The function uses the standard FFT definition and output values may grow by a
factor of <code>fftLen</code> when computing the forward transform. The
inverse transform includes a scale of <code>1/fftLen</code> as part of the
calculation and this matches the textbook definition of the inverse FFT.
@par
Pre-initialized data structures containing twiddle factors and bit reversal
tables are provided and defined in <code>arm_const_structs.h</code>. Include
this header in your function and then pass one of the constant structures as
an argument to arm_cfft_q31. For example:
@par
<code>arm_cfft_q31(arm_cfft_sR_q31_len64, pSrc, 1, 1)</code>
@par
computes a 64-point inverse complex FFT including bit reversal.
The data structures are treated as constant data and not modified during the
calculation. The same data structure can be reused for multiple transforms
including mixing forward and inverse transforms.
@par
Earlier releases of the library provided separate radix-2 and radix-4
algorithms that operated on floating-point data. These functions are still
provided but are deprecated. The older functions are slower and less general
than the new functions.
@par
An example of initialization of the constants for the arm_cfft_q31 function follows:
@code
const static arm_cfft_instance_q31 *S;
...
switch (length) {
case 16:
S = &arm_cfft_sR_q31_len16;
break;
case 32:
S = &arm_cfft_sR_q31_len32;
break;
case 64:
S = &arm_cfft_sR_q31_len64;
break;
case 128:
S = &arm_cfft_sR_q31_len128;
break;
case 256:
S = &arm_cfft_sR_q31_len256;
break;
case 512:
S = &arm_cfft_sR_q31_len512;
break;
case 1024:
S = &arm_cfft_sR_q31_len1024;
break;
case 2048:
S = &arm_cfft_sR_q31_len2048;
break;
case 4096:
S = &arm_cfft_sR_q31_len4096;
break;
}
@endcode
*/
void arm_cfft_radix8by2_f32 (arm_cfft_instance_f32 * S, float32_t * p1)
{
uint32_t L = S->fftLen;
float32_t * pCol1, * pCol2, * pMid1, * pMid2;
float32_t * p2 = p1 + L;
const float32_t * tw = (float32_t *) S->pTwiddle;
float32_t t1[4], t2[4], t3[4], t4[4], twR, twI;
float32_t m0, m1, m2, m3;
uint32_t l;
pCol1 = p1;
pCol2 = p2;
/* Define new length */
L >>= 1;
/* Initialize mid pointers */
pMid1 = p1 + L;
pMid2 = p2 + L;
/* do two dot Fourier transform */
for (l = L >> 2; l > 0; l-- )
{
t1[0] = p1[0];
t1[1] = p1[1];
t1[2] = p1[2];
t1[3] = p1[3];
t2[0] = p2[0];
t2[1] = p2[1];
t2[2] = p2[2];
t2[3] = p2[3];
t3[0] = pMid1[0];
t3[1] = pMid1[1];
t3[2] = pMid1[2];
t3[3] = pMid1[3];
t4[0] = pMid2[0];
t4[1] = pMid2[1];
t4[2] = pMid2[2];
t4[3] = pMid2[3];
*p1++ = t1[0] + t2[0];
*p1++ = t1[1] + t2[1];
*p1++ = t1[2] + t2[2];
*p1++ = t1[3] + t2[3]; /* col 1 */
t2[0] = t1[0] - t2[0];
t2[1] = t1[1] - t2[1];
t2[2] = t1[2] - t2[2];
t2[3] = t1[3] - t2[3]; /* for col 2 */
*pMid1++ = t3[0] + t4[0];
*pMid1++ = t3[1] + t4[1];
*pMid1++ = t3[2] + t4[2];
*pMid1++ = t3[3] + t4[3]; /* col 1 */
t4[0] = t4[0] - t3[0];
t4[1] = t4[1] - t3[1];
t4[2] = t4[2] - t3[2];
t4[3] = t4[3] - t3[3]; /* for col 2 */
twR = *tw++;
twI = *tw++;
/* multiply by twiddle factors */
m0 = t2[0] * twR;
m1 = t2[1] * twI;
m2 = t2[1] * twR;
m3 = t2[0] * twI;
/* R = R * Tr - I * Ti */
*p2++ = m0 + m1;
/* I = I * Tr + R * Ti */
*p2++ = m2 - m3;
/* use vertical symmetry */
/* 0.9988 - 0.0491i <==> -0.0491 - 0.9988i */
m0 = t4[0] * twI;
m1 = t4[1] * twR;
m2 = t4[1] * twI;
m3 = t4[0] * twR;
*pMid2++ = m0 - m1;
*pMid2++ = m2 + m3;
twR = *tw++;
twI = *tw++;
m0 = t2[2] * twR;
m1 = t2[3] * twI;
m2 = t2[3] * twR;
m3 = t2[2] * twI;
*p2++ = m0 + m1;
*p2++ = m2 - m3;
m0 = t4[2] * twI;
m1 = t4[3] * twR;
m2 = t4[3] * twI;
m3 = t4[2] * twR;
*pMid2++ = m0 - m1;
*pMid2++ = m2 + m3;
}
/* first col */
arm_radix8_butterfly_f32 (pCol1, L, (float32_t *) S->pTwiddle, 2U);
/* second col */
arm_radix8_butterfly_f32 (pCol2, L, (float32_t *) S->pTwiddle, 2U);
}
void arm_cfft_radix8by4_f32 (arm_cfft_instance_f32 * S, float32_t * p1)
{
uint32_t L = S->fftLen >> 1;
float32_t * pCol1, *pCol2, *pCol3, *pCol4, *pEnd1, *pEnd2, *pEnd3, *pEnd4;
const float32_t *tw2, *tw3, *tw4;
float32_t * p2 = p1 + L;
float32_t * p3 = p2 + L;
float32_t * p4 = p3 + L;
float32_t t2[4], t3[4], t4[4], twR, twI;
float32_t p1ap3_0, p1sp3_0, p1ap3_1, p1sp3_1;
float32_t m0, m1, m2, m3;
uint32_t l, twMod2, twMod3, twMod4;
pCol1 = p1; /* points to real values by default */
pCol2 = p2;
pCol3 = p3;
pCol4 = p4;
pEnd1 = p2 - 1; /* points to imaginary values by default */
pEnd2 = p3 - 1;
pEnd3 = p4 - 1;
pEnd4 = pEnd3 + L;
tw2 = tw3 = tw4 = (float32_t *) S->pTwiddle;
L >>= 1;
/* do four dot Fourier transform */
twMod2 = 2;
twMod3 = 4;
twMod4 = 6;
/* TOP */
p1ap3_0 = p1[0] + p3[0];
p1sp3_0 = p1[0] - p3[0];
p1ap3_1 = p1[1] + p3[1];
p1sp3_1 = p1[1] - p3[1];
/* col 2 */
t2[0] = p1sp3_0 + p2[1] - p4[1];
t2[1] = p1sp3_1 - p2[0] + p4[0];
/* col 3 */
t3[0] = p1ap3_0 - p2[0] - p4[0];
t3[1] = p1ap3_1 - p2[1] - p4[1];
/* col 4 */
t4[0] = p1sp3_0 - p2[1] + p4[1];
t4[1] = p1sp3_1 + p2[0] - p4[0];
/* col 1 */
*p1++ = p1ap3_0 + p2[0] + p4[0];
*p1++ = p1ap3_1 + p2[1] + p4[1];
/* Twiddle factors are ones */
*p2++ = t2[0];
*p2++ = t2[1];
*p3++ = t3[0];
*p3++ = t3[1];
*p4++ = t4[0];
*p4++ = t4[1];
tw2 += twMod2;
tw3 += twMod3;
tw4 += twMod4;
for (l = (L - 2) >> 1; l > 0; l-- )
{
/* TOP */
p1ap3_0 = p1[0] + p3[0];
p1sp3_0 = p1[0] - p3[0];
p1ap3_1 = p1[1] + p3[1];
p1sp3_1 = p1[1] - p3[1];
/* col 2 */
t2[0] = p1sp3_0 + p2[1] - p4[1];
t2[1] = p1sp3_1 - p2[0] + p4[0];
/* col 3 */
t3[0] = p1ap3_0 - p2[0] - p4[0];
t3[1] = p1ap3_1 - p2[1] - p4[1];
/* col 4 */
t4[0] = p1sp3_0 - p2[1] + p4[1];
t4[1] = p1sp3_1 + p2[0] - p4[0];
/* col 1 - top */
*p1++ = p1ap3_0 + p2[0] + p4[0];
*p1++ = p1ap3_1 + p2[1] + p4[1];
/* BOTTOM */
p1ap3_1 = pEnd1[-1] + pEnd3[-1];
p1sp3_1 = pEnd1[-1] - pEnd3[-1];
p1ap3_0 = pEnd1[ 0] + pEnd3[0];
p1sp3_0 = pEnd1[ 0] - pEnd3[0];
/* col 2 */
t2[2] = pEnd2[0] - pEnd4[0] + p1sp3_1;
t2[3] = pEnd1[0] - pEnd3[0] - pEnd2[-1] + pEnd4[-1];
/* col 3 */
t3[2] = p1ap3_1 - pEnd2[-1] - pEnd4[-1];
t3[3] = p1ap3_0 - pEnd2[ 0] - pEnd4[ 0];
/* col 4 */
t4[2] = pEnd2[ 0] - pEnd4[ 0] - p1sp3_1;
t4[3] = pEnd4[-1] - pEnd2[-1] - p1sp3_0;
/* col 1 - Bottom */
*pEnd1-- = p1ap3_0 + pEnd2[ 0] + pEnd4[ 0];
*pEnd1-- = p1ap3_1 + pEnd2[-1] + pEnd4[-1];
/* COL 2 */
/* read twiddle factors */
twR = *tw2++;
twI = *tw2++;
/* multiply by twiddle factors */
/* let Z1 = a + i(b), Z2 = c + i(d) */
/* => Z1 * Z2 = (a*c - b*d) + i(b*c + a*d) */
/* Top */
m0 = t2[0] * twR;
m1 = t2[1] * twI;
m2 = t2[1] * twR;
m3 = t2[0] * twI;
*p2++ = m0 + m1;
*p2++ = m2 - m3;
/* use vertical symmetry col 2 */
/* 0.9997 - 0.0245i <==> 0.0245 - 0.9997i */
/* Bottom */
m0 = t2[3] * twI;
m1 = t2[2] * twR;
m2 = t2[2] * twI;
m3 = t2[3] * twR;
*pEnd2-- = m0 - m1;
*pEnd2-- = m2 + m3;
/* COL 3 */
twR = tw3[0];
twI = tw3[1];
tw3 += twMod3;
/* Top */
m0 = t3[0] * twR;
m1 = t3[1] * twI;
m2 = t3[1] * twR;
m3 = t3[0] * twI;
*p3++ = m0 + m1;
*p3++ = m2 - m3;
/* use vertical symmetry col 3 */
/* 0.9988 - 0.0491i <==> -0.9988 - 0.0491i */
/* Bottom */
m0 = -t3[3] * twR;
m1 = t3[2] * twI;
m2 = t3[2] * twR;
m3 = t3[3] * twI;
*pEnd3-- = m0 - m1;
*pEnd3-- = m3 - m2;
/* COL 4 */
twR = tw4[0];
twI = tw4[1];
tw4 += twMod4;
/* Top */
m0 = t4[0] * twR;
m1 = t4[1] * twI;
m2 = t4[1] * twR;
m3 = t4[0] * twI;
*p4++ = m0 + m1;
*p4++ = m2 - m3;
/* use vertical symmetry col 4 */
/* 0.9973 - 0.0736i <==> -0.0736 + 0.9973i */
/* Bottom */
m0 = t4[3] * twI;
m1 = t4[2] * twR;
m2 = t4[2] * twI;
m3 = t4[3] * twR;
*pEnd4-- = m0 - m1;
*pEnd4-- = m2 + m3;
}
/* MIDDLE */
/* Twiddle factors are */
/* 1.0000 0.7071-0.7071i -1.0000i -0.7071-0.7071i */
p1ap3_0 = p1[0] + p3[0];
p1sp3_0 = p1[0] - p3[0];
p1ap3_1 = p1[1] + p3[1];
p1sp3_1 = p1[1] - p3[1];
/* col 2 */
t2[0] = p1sp3_0 + p2[1] - p4[1];
t2[1] = p1sp3_1 - p2[0] + p4[0];
/* col 3 */
t3[0] = p1ap3_0 - p2[0] - p4[0];
t3[1] = p1ap3_1 - p2[1] - p4[1];
/* col 4 */
t4[0] = p1sp3_0 - p2[1] + p4[1];
t4[1] = p1sp3_1 + p2[0] - p4[0];
/* col 1 - Top */
*p1++ = p1ap3_0 + p2[0] + p4[0];
*p1++ = p1ap3_1 + p2[1] + p4[1];
/* COL 2 */
twR = tw2[0];
twI = tw2[1];
m0 = t2[0] * twR;
m1 = t2[1] * twI;
m2 = t2[1] * twR;
m3 = t2[0] * twI;
*p2++ = m0 + m1;
*p2++ = m2 - m3;
/* COL 3 */
twR = tw3[0];
twI = tw3[1];
m0 = t3[0] * twR;
m1 = t3[1] * twI;
m2 = t3[1] * twR;
m3 = t3[0] * twI;
*p3++ = m0 + m1;
*p3++ = m2 - m3;
/* COL 4 */
twR = tw4[0];
twI = tw4[1];
m0 = t4[0] * twR;
m1 = t4[1] * twI;
m2 = t4[1] * twR;
m3 = t4[0] * twI;
*p4++ = m0 + m1;
*p4++ = m2 - m3;
/* first col */
arm_radix8_butterfly_f32 (pCol1, L, (float32_t *) S->pTwiddle, 4U);
/* second col */
arm_radix8_butterfly_f32 (pCol2, L, (float32_t *) S->pTwiddle, 4U);
/* third col */
arm_radix8_butterfly_f32 (pCol3, L, (float32_t *) S->pTwiddle, 4U);
/* fourth col */
arm_radix8_butterfly_f32 (pCol4, L, (float32_t *) S->pTwiddle, 4U);
}
/**
@addtogroup ComplexFFT
@{
*/
/**
@brief Processing function for the floating-point complex FFT.
@param[in] S points to an instance of the floating-point CFFT structure
@param[in,out] p1 points to the complex data buffer of size <code>2*fftLen</code>. Processing occurs in-place
@param[in] ifftFlag flag that selects transform direction
- value = 0: forward transform
- value = 1: inverse transform
@param[in] bitReverseFlag flag that enables / disables bit reversal of output
- value = 0: disables bit reversal of output
- value = 1: enables bit reversal of output
@return none
*/
void arm_cfft_f32(
const arm_cfft_instance_f32 * S,
float32_t * p1,
uint8_t ifftFlag,
uint8_t bitReverseFlag)
{
uint32_t L = S->fftLen, l;
float32_t invL, * pSrc;
if (ifftFlag == 1U)
{
/* Conjugate input data */
pSrc = p1 + 1;
for (l = 0; l < L; l++)
{
*pSrc = -*pSrc;
pSrc += 2;
}
}
switch (L)
{
case 16:
case 128:
case 1024:
arm_cfft_radix8by2_f32 ( (arm_cfft_instance_f32 *) S, p1);
break;
case 32:
case 256:
case 2048:
arm_cfft_radix8by4_f32 ( (arm_cfft_instance_f32 *) S, p1);
break;
case 64:
case 512:
case 4096:
arm_radix8_butterfly_f32 ( p1, L, (float32_t *) S->pTwiddle, 1);
break;
}
if ( bitReverseFlag )
arm_bitreversal_32 ((uint32_t*) p1, S->bitRevLength, S->pBitRevTable);
if (ifftFlag == 1U)
{
invL = 1.0f / (float32_t)L;
/* Conjugate and scale output data */
pSrc = p1;
for (l= 0; l < L; l++)
{
*pSrc++ *= invL ;
*pSrc = -(*pSrc) * invL;
pSrc++;
}
}
}
#endif /* defined(ARM_MATH_MVEF) && !defined(ARM_MATH_AUTOVECTORIZE) */
/**
@} end of ComplexFFT group
*/
Reference in New Issue View Git Blame Copy Permalink