arm_dct4_f32.c
1 /* ---------------------------------------------------------------------- 2 * Project: CMSIS DSP Library 3 * Title: arm_dct4_f32.c 4 * Description: Processing function of DCT4 & IDCT4 F32 5 * 6 * $Date: 23 April 2021 7 * $Revision: V1.9.0 8 * 9 * Target Processor: Cortex-M and Cortex-A cores 10 * -------------------------------------------------------------------- */ 11 /* 12 * Copyright (C) 2010-2021 ARM Limited or its affiliates. All rights reserved. 13 * 14 * SPDX-License-Identifier: Apache-2.0 15 * 16 * Licensed under the Apache License, Version 2.0 (the License); you may 17 * not use this file except in compliance with the License. 18 * You may obtain a copy of the License at 19 * 20 * www.apache.org/licenses/LICENSE-2.0 21 * 22 * Unless required by applicable law or agreed to in writing, software 23 * distributed under the License is distributed on an AS IS BASIS, WITHOUT 24 * WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. 25 * See the License for the specific language governing permissions and 26 * limitations under the License. 27 */ 28 29 #include "dsp/transform_functions.h" 30 31 /** 32 @ingroup groupTransforms 33 */ 34 35 /** 36 @defgroup DCT4_IDCT4 DCT Type IV Functions 37 38 Representation of signals by minimum number of values is important for storage and transmission. 39 The possibility of large discontinuity between the beginning and end of a period of a signal 40 in DFT can be avoided by extending the signal so that it is even-symmetric. 41 Discrete Cosine Transform (DCT) is constructed such that its energy is heavily concentrated in the lower part of the 42 spectrum and is very widely used in signal and image coding applications. 43 The family of DCTs (DCT type- 1,2,3,4) is the outcome of different combinations of homogeneous boundary conditions. 44 DCT has an excellent energy-packing capability, hence has many applications and in data compression in particular. 45 46 DCT is essentially the Discrete Fourier Transform(DFT) of an even-extended real signal. 47 Reordering of the input data makes the computation of DCT just a problem of 48 computing the DFT of a real signal with a few additional operations. 49 This approach provides regular, simple, and very efficient DCT algorithms for practical hardware and software implementations. 50 51 DCT type-II can be implemented using Fast fourier transform (FFT) internally, as the transform is applied on real values, Real FFT can be used. 52 DCT4 is implemented using DCT2 as their implementations are similar except with some added pre-processing and post-processing. 53 DCT2 implementation can be described in the following steps: 54 - Re-ordering input 55 - Calculating Real FFT 56 - Multiplication of weights and Real FFT output and getting real part from the product. 57 58 This process is explained by the block diagram below: 59 \image html DCT4.gif "Discrete Cosine Transform - type-IV" 60 61 @par Algorithm 62 The N-point type-IV DCT is defined as a real, linear transformation by the formula: 63 \image html DCT4Equation.gif 64 where <code>k = 0, 1, 2, ..., N-1</code> 65 @par 66 Its inverse is defined as follows: 67 \image html IDCT4Equation.gif 68 where <code>n = 0, 1, 2, ..., N-1</code> 69 @par 70 The DCT4 matrices become involutory (i.e. they are self-inverse) by multiplying with an overall scale factor of sqrt(2/N). 71 The symmetry of the transform matrix indicates that the fast algorithms for the forward 72 and inverse transform computation are identical. 73 Note that the implementation of Inverse DCT4 and DCT4 is same, hence same process function can be used for both. 74 75 @par Lengths supported by the transform: 76 As DCT4 internally uses Real FFT, it supports all the lengths 128, 512, 2048 and 8192. 77 The library provides separate functions for Q15, Q31, and floating-point data types. 78 79 @par Instance Structure 80 The instances for Real FFT and FFT, cosine values table and twiddle factor table are stored in an instance data structure. 81 A separate instance structure must be defined for each transform. 82 There are separate instance structure declarations for each of the 3 supported data types. 83 84 @par Initialization Functions 85 There is also an associated initialization function for each data type. 86 The initialization function performs the following operations: 87 - Sets the values of the internal structure fields. 88 - Initializes Real FFT as its process function is used internally in DCT4, by calling \ref arm_rfft_init_f32(). 89 @par 90 Use of the initialization function is optional. 91 However, if the initialization function is used, then the instance structure cannot be placed into a const data section. 92 To place an instance structure into a const data section, the instance structure must be manually initialized. 93 Manually initialize the instance structure as follows: 94 <pre> 95 arm_dct4_instance_f32 S = {N, Nby2, normalize, pTwiddle, pCosFactor, pRfft, pCfft}; 96 arm_dct4_instance_q31 S = {N, Nby2, normalize, pTwiddle, pCosFactor, pRfft, pCfft}; 97 arm_dct4_instance_q15 S = {N, Nby2, normalize, pTwiddle, pCosFactor, pRfft, pCfft}; 98 </pre> 99 where \c N is the length of the DCT4; \c Nby2 is half of the length of the DCT4; 100 \c normalize is normalizing factor used and is equal to <code>sqrt(2/N)</code>; 101 \c pTwiddle points to the twiddle factor table; 102 \c pCosFactor points to the cosFactor table; 103 \c pRfft points to the real FFT instance; 104 \c pCfft points to the complex FFT instance; 105 The CFFT and RFFT structures also needs to be initialized, refer to arm_cfft_radix4_f32() 106 and arm_rfft_f32() respectively for details regarding static initialization. 107 108 @par Fixed-Point Behavior 109 Care must be taken when using the fixed-point versions of the DCT4 transform functions. 110 In particular, the overflow and saturation behavior of the accumulator used in each function must be considered. 111 Refer to the function specific documentation below for usage guidelines. 112 */ 113 114 /** 115 @addtogroup DCT4_IDCT4 116 @{ 117 */ 118 119 /** 120 @brief Processing function for the floating-point DCT4/IDCT4. 121 @param[in] S points to an instance of the floating-point DCT4/IDCT4 structure 122 @param[in] pState points to state buffer 123 @param[in,out] pInlineBuffer points to the in-place input and output buffer 124 @return none 125 */ 126 127 void arm_dct4_f32( 128 const arm_dct4_instance_f32 * S, 129 float32_t * pState, 130 float32_t * pInlineBuffer) 131 { 132 const float32_t *weights = S->pTwiddle; /* Pointer to the Weights table */ 133 const float32_t *cosFact = S->pCosFactor; /* Pointer to the cos factors table */ 134 float32_t *pS1, *pS2, *pbuff; /* Temporary pointers for input buffer and pState buffer */ 135 float32_t in; /* Temporary variable */ 136 uint32_t i; /* Loop counter */ 137 138 139 /* DCT4 computation involves DCT2 (which is calculated using RFFT) 140 * along with some pre-processing and post-processing. 141 * Computational procedure is explained as follows: 142 * (a) Pre-processing involves multiplying input with cos factor, 143 * r(n) = 2 * u(n) * cos(pi*(2*n+1)/(4*n)) 144 * where, 145 * r(n) -- output of preprocessing 146 * u(n) -- input to preprocessing(actual Source buffer) 147 * (b) Calculation of DCT2 using FFT is divided into three steps: 148 * Step1: Re-ordering of even and odd elements of input. 149 * Step2: Calculating FFT of the re-ordered input. 150 * Step3: Taking the real part of the product of FFT output and weights. 151 * (c) Post-processing - DCT4 can be obtained from DCT2 output using the following equation: 152 * Y4(k) = Y2(k) - Y4(k-1) and Y4(-1) = Y4(0) 153 * where, 154 * Y4 -- DCT4 output, Y2 -- DCT2 output 155 * (d) Multiplying the output with the normalizing factor sqrt(2/N). 156 */ 157 158 /*-------- Pre-processing ------------*/ 159 /* Multiplying input with cos factor i.e. r(n) = 2 * x(n) * cos(pi*(2*n+1)/(4*n)) */ 160 arm_scale_f32(pInlineBuffer, 2.0f, pInlineBuffer, S->N); 161 arm_mult_f32(pInlineBuffer, cosFact, pInlineBuffer, S->N); 162 163 /* ---------------------------------------------------------------- 164 * Step1: Re-ordering of even and odd elements as 165 * pState[i] = pInlineBuffer[2*i] and 166 * pState[N-i-1] = pInlineBuffer[2*i+1] where i = 0 to N/2 167 ---------------------------------------------------------------------*/ 168 169 /* pS1 initialized to pState */ 170 pS1 = pState; 171 172 /* pS2 initialized to pState+N-1, so that it points to the end of the state buffer */ 173 pS2 = pState + (S->N - 1U); 174 175 /* pbuff initialized to input buffer */ 176 pbuff = pInlineBuffer; 177 178 179 #if defined (ARM_MATH_LOOPUNROLL) 180 181 /* Initializing the loop counter to N/2 >> 2 for loop unrolling by 4 */ 182 i = S->Nby2 >> 2U; 183 184 /* First part of the processing with loop unrolling. Compute 4 outputs at a time. 185 ** a second loop below computes the remaining 1 to 3 samples. */ 186 do 187 { 188 /* Re-ordering of even and odd elements */ 189 /* pState[i] = pInlineBuffer[2*i] */ 190 *pS1++ = *pbuff++; 191 /* pState[N-i-1] = pInlineBuffer[2*i+1] */ 192 *pS2-- = *pbuff++; 193 194 *pS1++ = *pbuff++; 195 *pS2-- = *pbuff++; 196 197 *pS1++ = *pbuff++; 198 *pS2-- = *pbuff++; 199 200 *pS1++ = *pbuff++; 201 *pS2-- = *pbuff++; 202 203 /* Decrement loop counter */ 204 i--; 205 } while (i > 0U); 206 207 /* pbuff initialized to input buffer */ 208 pbuff = pInlineBuffer; 209 210 /* pS1 initialized to pState */ 211 pS1 = pState; 212 213 /* Initializing the loop counter to N/4 instead of N for loop unrolling */ 214 i = S->N >> 2U; 215 216 /* Processing with loop unrolling 4 times as N is always multiple of 4. 217 * Compute 4 outputs at a time */ 218 do 219 { 220 /* Writing the re-ordered output back to inplace input buffer */ 221 *pbuff++ = *pS1++; 222 *pbuff++ = *pS1++; 223 *pbuff++ = *pS1++; 224 *pbuff++ = *pS1++; 225 226 /* Decrement the loop counter */ 227 i--; 228 } while (i > 0U); 229 230 231 /* --------------------------------------------------------- 232 * Step2: Calculate RFFT for N-point input 233 * ---------------------------------------------------------- */ 234 /* pInlineBuffer is real input of length N , pState is the complex output of length 2N */ 235 arm_rfft_f32 (S->pRfft, pInlineBuffer, pState); 236 237 /*---------------------------------------------------------------------- 238 * Step3: Multiply the FFT output with the weights. 239 *----------------------------------------------------------------------*/ 240 arm_cmplx_mult_cmplx_f32 (pState, weights, pState, S->N); 241 242 /* ----------- Post-processing ---------- */ 243 /* DCT-IV can be obtained from DCT-II by the equation, 244 * Y4(k) = Y2(k) - Y4(k-1) and Y4(-1) = Y4(0) 245 * Hence, Y4(0) = Y2(0)/2 */ 246 /* Getting only real part from the output and Converting to DCT-IV */ 247 248 /* Initializing the loop counter to N >> 2 for loop unrolling by 4 */ 249 i = (S->N - 1U) >> 2U; 250 251 /* pbuff initialized to input buffer. */ 252 pbuff = pInlineBuffer; 253 254 /* pS1 initialized to pState */ 255 pS1 = pState; 256 257 /* Calculating Y4(0) from Y2(0) using Y4(0) = Y2(0)/2 */ 258 in = *pS1++ * (float32_t) 0.5; 259 /* input buffer acts as inplace, so output values are stored in the input itself. */ 260 *pbuff++ = in; 261 262 /* pState pointer is incremented twice as the real values are located alternatively in the array */ 263 pS1++; 264 265 /* First part of the processing with loop unrolling. Compute 4 outputs at a time. 266 ** a second loop below computes the remaining 1 to 3 samples. */ 267 do 268 { 269 /* Calculating Y4(1) to Y4(N-1) from Y2 using equation Y4(k) = Y2(k) - Y4(k-1) */ 270 /* pState pointer (pS1) is incremented twice as the real values are located alternatively in the array */ 271 in = *pS1++ - in; 272 *pbuff++ = in; 273 /* points to the next real value */ 274 pS1++; 275 276 in = *pS1++ - in; 277 *pbuff++ = in; 278 pS1++; 279 280 in = *pS1++ - in; 281 *pbuff++ = in; 282 pS1++; 283 284 in = *pS1++ - in; 285 *pbuff++ = in; 286 pS1++; 287 288 /* Decrement the loop counter */ 289 i--; 290 } while (i > 0U); 291 292 /* If the blockSize is not a multiple of 4, compute any remaining output samples here. 293 ** No loop unrolling is used. */ 294 i = (S->N - 1U) % 0x4U; 295 296 while (i > 0U) 297 { 298 /* Calculating Y4(1) to Y4(N-1) from Y2 using equation Y4(k) = Y2(k) - Y4(k-1) */ 299 /* pState pointer (pS1) is incremented twice as the real values are located alternatively in the array */ 300 in = *pS1++ - in; 301 *pbuff++ = in; 302 303 /* points to the next real value */ 304 pS1++; 305 306 /* Decrement the loop counter */ 307 i--; 308 } 309 310 311 /*------------ Normalizing the output by multiplying with the normalizing factor ----------*/ 312 313 /* Initializing the loop counter to N/4 instead of N for loop unrolling */ 314 i = S->N >> 2U; 315 316 /* pbuff initialized to the pInlineBuffer(now contains the output values) */ 317 pbuff = pInlineBuffer; 318 319 /* Processing with loop unrolling 4 times as N is always multiple of 4. Compute 4 outputs at a time */ 320 do 321 { 322 /* Multiplying pInlineBuffer with the normalizing factor sqrt(2/N) */ 323 in = *pbuff; 324 *pbuff++ = in * S->normalize; 325 326 in = *pbuff; 327 *pbuff++ = in * S->normalize; 328 329 in = *pbuff; 330 *pbuff++ = in * S->normalize; 331 332 in = *pbuff; 333 *pbuff++ = in * S->normalize; 334 335 /* Decrement the loop counter */ 336 i--; 337 } while (i > 0U); 338 339 340 #else 341 342 /* Initializing the loop counter to N/2 */ 343 i = S->Nby2; 344 345 do 346 { 347 /* Re-ordering of even and odd elements */ 348 /* pState[i] = pInlineBuffer[2*i] */ 349 *pS1++ = *pbuff++; 350 /* pState[N-i-1] = pInlineBuffer[2*i+1] */ 351 *pS2-- = *pbuff++; 352 353 /* Decrement the loop counter */ 354 i--; 355 } while (i > 0U); 356 357 /* pbuff initialized to input buffer */ 358 pbuff = pInlineBuffer; 359 360 /* pS1 initialized to pState */ 361 pS1 = pState; 362 363 /* Initializing the loop counter */ 364 i = S->N; 365 366 do 367 { 368 /* Writing the re-ordered output back to inplace input buffer */ 369 *pbuff++ = *pS1++; 370 371 /* Decrement the loop counter */ 372 i--; 373 } while (i > 0U); 374 375 376 /* --------------------------------------------------------- 377 * Step2: Calculate RFFT for N-point input 378 * ---------------------------------------------------------- */ 379 /* pInlineBuffer is real input of length N , pState is the complex output of length 2N */ 380 arm_rfft_f32 (S->pRfft, pInlineBuffer, pState); 381 382 /*---------------------------------------------------------------------- 383 * Step3: Multiply the FFT output with the weights. 384 *----------------------------------------------------------------------*/ 385 arm_cmplx_mult_cmplx_f32 (pState, weights, pState, S->N); 386 387 /* ----------- Post-processing ---------- */ 388 /* DCT-IV can be obtained from DCT-II by the equation, 389 * Y4(k) = Y2(k) - Y4(k-1) and Y4(-1) = Y4(0) 390 * Hence, Y4(0) = Y2(0)/2 */ 391 /* Getting only real part from the output and Converting to DCT-IV */ 392 393 /* pbuff initialized to input buffer. */ 394 pbuff = pInlineBuffer; 395 396 /* pS1 initialized to pState */ 397 pS1 = pState; 398 399 /* Calculating Y4(0) from Y2(0) using Y4(0) = Y2(0)/2 */ 400 in = *pS1++ * (float32_t) 0.5; 401 /* input buffer acts as inplace, so output values are stored in the input itself. */ 402 *pbuff++ = in; 403 404 /* pState pointer is incremented twice as the real values are located alternatively in the array */ 405 pS1++; 406 407 /* Initializing the loop counter */ 408 i = (S->N - 1U); 409 410 do 411 { 412 /* Calculating Y4(1) to Y4(N-1) from Y2 using equation Y4(k) = Y2(k) - Y4(k-1) */ 413 /* pState pointer (pS1) is incremented twice as the real values are located alternatively in the array */ 414 in = *pS1++ - in; 415 *pbuff++ = in; 416 417 /* points to the next real value */ 418 pS1++; 419 420 /* Decrement loop counter */ 421 i--; 422 } while (i > 0U); 423 424 /*------------ Normalizing the output by multiplying with the normalizing factor ----------*/ 425 426 /* Initializing loop counter */ 427 i = S->N; 428 429 /* pbuff initialized to the pInlineBuffer (now contains the output values) */ 430 pbuff = pInlineBuffer; 431 432 do 433 { 434 /* Multiplying pInlineBuffer with the normalizing factor sqrt(2/N) */ 435 in = *pbuff; 436 *pbuff++ = in * S->normalize; 437 438 /* Decrement loop counter */ 439 i--; 440 } while (i > 0U); 441 442 #endif /* #if defined (ARM_MATH_LOOPUNROLL) */ 443 444 } 445 446 /** 447 @} end of DCT4_IDCT4 group 448 */