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
arm_dct4_f32(const arm_dct4_instance_f32 * S,float32_t * pState,float32_t * pInlineBuffer)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 */
449
