/* ----------------------------------------------------------------------
* Project: CMSIS DSP Library
* Title: arm_rfft_f32.c
* Description: RFFT & RIFFT Floating point process function
*
* $Date: 18. March 2019
* $Revision: V1.6.0
*
* Target Processor: Cortex-M cores
* -------------------------------------------------------------------- */
/*
* Copyright (C) 2010-2019 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 "arm_math.h"
void stage_rfft_f32(
const arm_rfft_fast_instance_f32 * S,
float32_t * p,
float32_t * pOut)
{
uint32_t k; /* Loop Counter */
float32_t twR, twI; /* RFFT Twiddle coefficients */
const float32_t * pCoeff = S->pTwiddleRFFT; /* Points to RFFT Twiddle factors */
float32_t *pA = p; /* increasing pointer */
float32_t *pB = p; /* decreasing pointer */
float32_t xAR, xAI, xBR, xBI; /* temporary variables */
float32_t t1a, t1b; /* temporary variables */
float32_t p0, p1, p2, p3; /* temporary variables */
k = (S->Sint).fftLen - 1;
/* Pack first and last sample of the frequency domain together */
xBR = pB[0];
xBI = pB[1];
xAR = pA[0];
xAI = pA[1];
twR = *pCoeff++ ;
twI = *pCoeff++ ;
// U1 = XA(1) + XB(1); % It is real
t1a = xBR + xAR ;
// U2 = XB(1) - XA(1); % It is imaginary
t1b = xBI + xAI ;
// real(tw * (xB - xA)) = twR * (xBR - xAR) - twI * (xBI - xAI);
// imag(tw * (xB - xA)) = twI * (xBR - xAR) + twR * (xBI - xAI);
*pOut++ = 0.5f * ( t1a + t1b );
*pOut++ = 0.5f * ( t1a - t1b );
// XA(1) = 1/2*( U1 - imag(U2) + i*( U1 +imag(U2) ));
pB = p + 2*k;
pA += 2;
do
{
/*
function X = my_split_rfft(X, ifftFlag)
% X is a series of real numbers
L = length(X);
XC = X(1:2:end) +i*X(2:2:end);
XA = fft(XC);
XB = conj(XA([1 end:-1:2]));
TW = i*exp(-2*pi*i*[0:L/2-1]/L).';
for l = 2:L/2
XA(l) = 1/2 * (XA(l) + XB(l) + TW(l) * (XB(l) - XA(l)));
end
XA(1) = 1/2* (XA(1) + XB(1) + TW(1) * (XB(1) - XA(1))) + i*( 1/2*( XA(1) + XB(1) + i*( XA(1) - XB(1))));
X = XA;
*/
xBI = pB[1];
xBR = pB[0];
xAR = pA[0];
xAI = pA[1];
twR = *pCoeff++;
twI = *pCoeff++;
t1a = xBR - xAR ;
t1b = xBI + xAI ;
// real(tw * (xB - xA)) = twR * (xBR - xAR) - twI * (xBI - xAI);
// imag(tw * (xB - xA)) = twI * (xBR - xAR) + twR * (xBI - xAI);
p0 = twR * t1a;
p1 = twI * t1a;
p2 = twR * t1b;
p3 = twI * t1b;
*pOut++ = 0.5f * (xAR + xBR + p0 + p3 ); //xAR
*pOut++ = 0.5f * (xAI - xBI + p1 - p2 ); //xAI
pA += 2;
pB -= 2;
k--;
} while (k > 0U);
}
/* Prepares data for inverse cfft */
void merge_rfft_f32(
const arm_rfft_fast_instance_f32 * S,
float32_t * p,
float32_t * pOut)
{
uint32_t k; /* Loop Counter */
float32_t twR, twI; /* RFFT Twiddle coefficients */
const float32_t *pCoeff = S->pTwiddleRFFT; /* Points to RFFT Twiddle factors */
float32_t *pA = p; /* increasing pointer */
float32_t *pB = p; /* decreasing pointer */
float32_t xAR, xAI, xBR, xBI; /* temporary variables */
float32_t t1a, t1b, r, s, t, u; /* temporary variables */
k = (S->Sint).fftLen - 1;
xAR = pA[0];
xAI = pA[1];
pCoeff += 2 ;
*pOut++ = 0.5f * ( xAR + xAI );
*pOut++ = 0.5f * ( xAR - xAI );
pB = p + 2*k ;
pA += 2 ;
while (k > 0U)
{
/* G is half of the frequency complex spectrum */
//for k = 2:N
// Xk(k) = 1/2 * (G(k) + conj(G(N-k+2)) + Tw(k)*( G(k) - conj(G(N-k+2))));
xBI = pB[1] ;
xBR = pB[0] ;
xAR = pA[0];
xAI = pA[1];
twR = *pCoeff++;
twI = *pCoeff++;
t1a = xAR - xBR ;
t1b = xAI + xBI ;
r = twR * t1a;
s = twI * t1b;
t = twI * t1a;
u = twR * t1b;
// real(tw * (xA - xB)) = twR * (xAR - xBR) - twI * (xAI - xBI);
// imag(tw * (xA - xB)) = twI * (xAR - xBR) + twR * (xAI - xBI);
*pOut++ = 0.5f * (xAR + xBR - r - s ); //xAR
*pOut++ = 0.5f * (xAI - xBI + t - u ); //xAI
pA += 2;
pB -= 2;
k--;
}
}
/**
@ingroup groupTransforms
*/
/**
@defgroup RealFFT Real FFT Functions
@par
The CMSIS DSP library includes specialized algorithms for computing the
FFT of real data sequences. The FFT is defined over complex data but
in many applications the input is real. Real FFT algorithms take advantage
of the symmetry properties of the FFT and have a speed advantage over complex
algorithms of the same length.
@par
The Fast RFFT algorith relays on the mixed radix CFFT that save processor usage.
@par
The real length N forward FFT of a sequence is computed using the steps shown below.
@par
\image html RFFT.gif "Real Fast Fourier Transform"
@par
The real sequence is initially treated as if it were complex to perform a CFFT.
Later, a processing stage reshapes the data to obtain half of the frequency spectrum
in complex format. Except the first complex number that contains the two real numbers
X[0] and X[N/2] all the data is complex. In other words, the first complex sample
contains two real values packed.
@par
The input for the inverse RFFT should keep the same format as the output of the
forward RFFT. A first processing stage pre-process the data to later perform an
inverse CFFT.
@par
\image html RIFFT.gif "Real Inverse Fast Fourier Transform"
@par
The algorithms for floating-point, Q15, and Q31 data are slightly different
and we describe each algorithm in turn.
@par Floating-point
The main functions are \ref arm_rfft_fast_f32() and \ref arm_rfft_fast_init_f32().
The older functions \ref arm_rfft_f32() and \ref arm_rfft_init_f32() have been deprecated
but are still documented.
@par
The FFT of a real N-point sequence has even symmetry in the frequency domain.
The second half of the data equals the conjugate of the first half flipped in frequency.
Looking at the data, we see that we can uniquely represent the FFT using only N/2 complex numbers.
These are packed into the output array in alternating real and imaginary components:
@par
X = { real[0], imag[0], real[1], imag[1], real[2], imag[2] ...
real[(N/2)-1], imag[(N/2)-1 }
@par
It happens that the first complex number (real[0], imag[0]) is actually
all real. real[0] represents the DC offset, and imag[0] should be 0.
(real[1], imag[1]) is the fundamental frequency, (real[2], imag[2]) is
the first harmonic and so on.
@par
The real FFT functions pack the frequency domain data in this fashion.
The forward transform outputs the data in this form and the inverse
transform expects input data in this form. The function always performs
the needed bitreversal so that the input and output data is always in
normal order. The functions support lengths of [32, 64, 128, ..., 4096]
samples.
@par Q15 and Q31
The real algorithms are defined in a similar manner and utilize N/2 complex
transforms behind the scenes.
@par
The complex transforms used internally include scaling to prevent fixed-point
overflows. The overall scaling equals 1/(fftLen/2).
@par
A separate instance structure must be defined for each transform used but
twiddle factor and bit reversal tables can be reused.
@par
There is also an associated initialization function for each data type.
The initialization function performs the following operations:
- Sets the values of the internal structure fields.
- Initializes twiddle factor table and bit reversal table pointers.
- Initializes the internal complex FFT data structure.
@par
Use of the initialization function is optional.
However, if the initialization function is used, then the instance structure
cannot be placed into a const data section. To place an instance structure
into a const data section, the instance structure should be manually
initialized as follows:
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};
where fftLenReal is the length of the real transform;
fftLenBy2 length of the internal complex transform.
ifftFlagR Selects forward (=0) or inverse (=1) transform.
bitReverseFlagR Selects bit reversed output (=0) or normal order
output (=1).
twidCoefRModifier stride modifier for the twiddle factor table.
The value is based on the FFT length;
pTwiddleARealpoints to the A array of twiddle coefficients;
pTwiddleBRealpoints to the B array of twiddle coefficients;
pCfft points to the CFFT Instance structure. The CFFT structure
must also be initialized. Refer to arm_cfft_radix4_f32() for details regarding
static initialization of the complex FFT instance structure.
*/
/**
@addtogroup RealFFT
@{
*/
/**
@brief Processing function for the floating-point real FFT.
@param[in] S points to an arm_rfft_fast_instance_f32 structure
@param[in] p points to input buffer
@param[in] pOut points to output buffer
@param[in] ifftFlag
- value = 0: RFFT
- value = 1: RIFFT
@return none
*/
void arm_rfft_fast_f32(
arm_rfft_fast_instance_f32 * S,
float32_t * p,
float32_t * pOut,
uint8_t ifftFlag)
{
arm_cfft_instance_f32 * Sint = &(S->Sint);
Sint->fftLen = S->fftLenRFFT / 2;
/* Calculation of Real FFT */
if (ifftFlag)
{
/* Real FFT compression */
merge_rfft_f32(S, p, pOut);
/* Complex radix-4 IFFT process */
arm_cfft_f32( Sint, pOut, ifftFlag, 1);
}
else
{
/* Calculation of RFFT of input */
arm_cfft_f32( Sint, p, ifftFlag, 1);
/* Real FFT extraction */
stage_rfft_f32(S, p, pOut);
}
}
/**
* @} end of RealFFT group
*/