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 \f[
64 X_c(k) = \sqrt{\frac{2}{N}}\sum_{n=0}^{N-1} x(n)cos\Big[\Big(n+\frac{1}{2}\Big)\Big(k+\frac{1}{2}\Big)\frac{\pi}{N}\Big]
65 \f]
66 where <code>k = 0, 1, 2, ..., N-1</code>
67 @par
68 Its inverse is defined as follows:
69 \f[
70 x(n) = \sqrt{\frac{2}{N}}\sum_{k=0}^{N-1} X_c(k)cos\Big[\Big(n+\frac{1}{2}\Big)\Big(k+\frac{1}{2}\Big)\frac{\pi}{N}\Big]
71 \f]
72 where <code>n = 0, 1, 2, ..., N-1</code>
73 @par
74 The DCT4 matrices become involutory (i.e. they are self-inverse) by multiplying with an overall scale factor of sqrt(2/N).
75 The symmetry of the transform matrix indicates that the fast algorithms for the forward
76 and inverse transform computation are identical.
77 Note that the implementation of Inverse DCT4 and DCT4 is same, hence same process function can be used for both.
78
79 @par Lengths supported by the transform:
80 As DCT4 internally uses Real FFT, it supports all the lengths 128, 512, 2048 and 8192.
81 The library provides separate functions for Q15, Q31, and floating-point data types.
82
83 @par Instance Structure
84 The instances for Real FFT and FFT, cosine values table and twiddle factor table are stored in an instance data structure.
85 A separate instance structure must be defined for each transform.
86 There are separate instance structure declarations for each of the 3 supported data types.
87
88 @par Initialization Functions
89 There is also an associated initialization function for each data type.
90 The initialization function performs the following operations:
91 - Sets the values of the internal structure fields.
92 - Initializes Real FFT as its process function is used internally in DCT4, by calling \ref arm_rfft_init_f32().
93 @par
94 Use of the initialization function is optional.
95 However, if the initialization function is used, then the instance structure cannot be placed into a const data section.
96 To place an instance structure into a const data section, the instance structure must be manually initialized.
97 Manually initialize the instance structure as follows:
98 <pre>
99 arm_dct4_instance_f32 S = {N, Nby2, normalize, pTwiddle, pCosFactor, pRfft, pCfft};
100 arm_dct4_instance_q31 S = {N, Nby2, normalize, pTwiddle, pCosFactor, pRfft, pCfft};
101 arm_dct4_instance_q15 S = {N, Nby2, normalize, pTwiddle, pCosFactor, pRfft, pCfft};
102 </pre>
103 where \c N is the length of the DCT4; \c Nby2 is half of the length of the DCT4;
104 \c normalize is normalizing factor used and is equal to <code>sqrt(2/N)</code>;
105 \c pTwiddle points to the twiddle factor table;
106 \c pCosFactor points to the cosFactor table;
107 \c pRfft points to the real FFT instance;
108 \c pCfft points to the complex FFT instance;
109 The CFFT and RFFT structures also needs to be initialized, refer to arm_cfft_radix4_f32()
110 and arm_rfft_f32() respectively for details regarding static initialization.
111
112 @par Fixed-Point Behavior
113 Care must be taken when using the fixed-point versions of the DCT4 transform functions.
114 In particular, the overflow and saturation behavior of the accumulator used in each function must be considered.
115 Refer to the function specific documentation below for usage guidelines.
116 */
117
118 /**
119 @addtogroup DCT4F32
120 @{
121 */
122
123 /**
124 @brief Processing function for the floating-point DCT4/IDCT4.
125 @deprecated Do not use this function. It is using a deprecated version of the RFFT.
126 @param[in] S points to an instance of the floating-point DCT4/IDCT4 structure
127 @param[in] pState points to state buffer
128 @param[in,out] pInlineBuffer points to the in-place input and output buffer
129 */
130
arm_dct4_f32(const arm_dct4_instance_f32 * S,float32_t * pState,float32_t * pInlineBuffer)131 void arm_dct4_f32(
132 const arm_dct4_instance_f32 * S,
133 float32_t * pState,
134 float32_t * pInlineBuffer)
135 {
136 const float32_t *weights = S->pTwiddle; /* Pointer to the Weights table */
137 const float32_t *cosFact = S->pCosFactor; /* Pointer to the cos factors table */
138 float32_t *pS1, *pS2, *pbuff; /* Temporary pointers for input buffer and pState buffer */
139 float32_t in; /* Temporary variable */
140 uint32_t i; /* Loop counter */
141
142
143 /* DCT4 computation involves DCT2 (which is calculated using RFFT)
144 * along with some pre-processing and post-processing.
145 * Computational procedure is explained as follows:
146 * (a) Pre-processing involves multiplying input with cos factor,
147 * r(n) = 2 * u(n) * cos(pi*(2*n+1)/(4*n))
148 * where,
149 * r(n) -- output of preprocessing
150 * u(n) -- input to preprocessing(actual Source buffer)
151 * (b) Calculation of DCT2 using FFT is divided into three steps:
152 * Step1: Re-ordering of even and odd elements of input.
153 * Step2: Calculating FFT of the re-ordered input.
154 * Step3: Taking the real part of the product of FFT output and weights.
155 * (c) Post-processing - DCT4 can be obtained from DCT2 output using the following equation:
156 * Y4(k) = Y2(k) - Y4(k-1) and Y4(-1) = Y4(0)
157 * where,
158 * Y4 -- DCT4 output, Y2 -- DCT2 output
159 * (d) Multiplying the output with the normalizing factor sqrt(2/N).
160 */
161
162 /*-------- Pre-processing ------------*/
163 /* Multiplying input with cos factor i.e. r(n) = 2 * x(n) * cos(pi*(2*n+1)/(4*n)) */
164 arm_scale_f32(pInlineBuffer, 2.0f, pInlineBuffer, S->N);
165 arm_mult_f32(pInlineBuffer, cosFact, pInlineBuffer, S->N);
166
167 /* ----------------------------------------------------------------
168 * Step1: Re-ordering of even and odd elements as
169 * pState[i] = pInlineBuffer[2*i] and
170 * pState[N-i-1] = pInlineBuffer[2*i+1] where i = 0 to N/2
171 ---------------------------------------------------------------------*/
172
173 /* pS1 initialized to pState */
174 pS1 = pState;
175
176 /* pS2 initialized to pState+N-1, so that it points to the end of the state buffer */
177 pS2 = pState + (S->N - 1U);
178
179 /* pbuff initialized to input buffer */
180 pbuff = pInlineBuffer;
181
182
183 #if defined (ARM_MATH_LOOPUNROLL)
184
185 /* Initializing the loop counter to N/2 >> 2 for loop unrolling by 4 */
186 i = S->Nby2 >> 2U;
187
188 /* First part of the processing with loop unrolling. Compute 4 outputs at a time.
189 ** a second loop below computes the remaining 1 to 3 samples. */
190 do
191 {
192 /* Re-ordering of even and odd elements */
193 /* pState[i] = pInlineBuffer[2*i] */
194 *pS1++ = *pbuff++;
195 /* pState[N-i-1] = pInlineBuffer[2*i+1] */
196 *pS2-- = *pbuff++;
197
198 *pS1++ = *pbuff++;
199 *pS2-- = *pbuff++;
200
201 *pS1++ = *pbuff++;
202 *pS2-- = *pbuff++;
203
204 *pS1++ = *pbuff++;
205 *pS2-- = *pbuff++;
206
207 /* Decrement loop counter */
208 i--;
209 } while (i > 0U);
210
211 /* pbuff initialized to input buffer */
212 pbuff = pInlineBuffer;
213
214 /* pS1 initialized to pState */
215 pS1 = pState;
216
217 /* Initializing the loop counter to N/4 instead of N for loop unrolling */
218 i = S->N >> 2U;
219
220 /* Processing with loop unrolling 4 times as N is always multiple of 4.
221 * Compute 4 outputs at a time */
222 do
223 {
224 /* Writing the re-ordered output back to inplace input buffer */
225 *pbuff++ = *pS1++;
226 *pbuff++ = *pS1++;
227 *pbuff++ = *pS1++;
228 *pbuff++ = *pS1++;
229
230 /* Decrement the loop counter */
231 i--;
232 } while (i > 0U);
233
234
235 /* ---------------------------------------------------------
236 * Step2: Calculate RFFT for N-point input
237 * ---------------------------------------------------------- */
238 /* pInlineBuffer is real input of length N , pState is the complex output of length 2N */
239 arm_rfft_f32 (S->pRfft, pInlineBuffer, pState);
240
241 /*----------------------------------------------------------------------
242 * Step3: Multiply the FFT output with the weights.
243 *----------------------------------------------------------------------*/
244 arm_cmplx_mult_cmplx_f32 (pState, weights, pState, S->N);
245
246 /* ----------- Post-processing ---------- */
247 /* DCT-IV can be obtained from DCT-II by the equation,
248 * Y4(k) = Y2(k) - Y4(k-1) and Y4(-1) = Y4(0)
249 * Hence, Y4(0) = Y2(0)/2 */
250 /* Getting only real part from the output and Converting to DCT-IV */
251
252 /* Initializing the loop counter to N >> 2 for loop unrolling by 4 */
253 i = (S->N - 1U) >> 2U;
254
255 /* pbuff initialized to input buffer. */
256 pbuff = pInlineBuffer;
257
258 /* pS1 initialized to pState */
259 pS1 = pState;
260
261 /* Calculating Y4(0) from Y2(0) using Y4(0) = Y2(0)/2 */
262 in = *pS1++ * (float32_t) 0.5;
263 /* input buffer acts as inplace, so output values are stored in the input itself. */
264 *pbuff++ = in;
265
266 /* pState pointer is incremented twice as the real values are located alternatively in the array */
267 pS1++;
268
269 /* First part of the processing with loop unrolling. Compute 4 outputs at a time.
270 ** a second loop below computes the remaining 1 to 3 samples. */
271 do
272 {
273 /* Calculating Y4(1) to Y4(N-1) from Y2 using equation Y4(k) = Y2(k) - Y4(k-1) */
274 /* pState pointer (pS1) is incremented twice as the real values are located alternatively in the array */
275 in = *pS1++ - in;
276 *pbuff++ = in;
277 /* points to the next real value */
278 pS1++;
279
280 in = *pS1++ - in;
281 *pbuff++ = in;
282 pS1++;
283
284 in = *pS1++ - in;
285 *pbuff++ = in;
286 pS1++;
287
288 in = *pS1++ - in;
289 *pbuff++ = in;
290 pS1++;
291
292 /* Decrement the loop counter */
293 i--;
294 } while (i > 0U);
295
296 /* If the blockSize is not a multiple of 4, compute any remaining output samples here.
297 ** No loop unrolling is used. */
298 i = (S->N - 1U) % 0x4U;
299
300 while (i > 0U)
301 {
302 /* Calculating Y4(1) to Y4(N-1) from Y2 using equation Y4(k) = Y2(k) - Y4(k-1) */
303 /* pState pointer (pS1) is incremented twice as the real values are located alternatively in the array */
304 in = *pS1++ - in;
305 *pbuff++ = in;
306
307 /* points to the next real value */
308 pS1++;
309
310 /* Decrement the loop counter */
311 i--;
312 }
313
314
315 /*------------ Normalizing the output by multiplying with the normalizing factor ----------*/
316
317 /* Initializing the loop counter to N/4 instead of N for loop unrolling */
318 i = S->N >> 2U;
319
320 /* pbuff initialized to the pInlineBuffer(now contains the output values) */
321 pbuff = pInlineBuffer;
322
323 /* Processing with loop unrolling 4 times as N is always multiple of 4. Compute 4 outputs at a time */
324 do
325 {
326 /* Multiplying pInlineBuffer with the normalizing factor sqrt(2/N) */
327 in = *pbuff;
328 *pbuff++ = in * S->normalize;
329
330 in = *pbuff;
331 *pbuff++ = in * S->normalize;
332
333 in = *pbuff;
334 *pbuff++ = in * S->normalize;
335
336 in = *pbuff;
337 *pbuff++ = in * S->normalize;
338
339 /* Decrement the loop counter */
340 i--;
341 } while (i > 0U);
342
343
344 #else
345
346 /* Initializing the loop counter to N/2 */
347 i = S->Nby2;
348
349 do
350 {
351 /* Re-ordering of even and odd elements */
352 /* pState[i] = pInlineBuffer[2*i] */
353 *pS1++ = *pbuff++;
354 /* pState[N-i-1] = pInlineBuffer[2*i+1] */
355 *pS2-- = *pbuff++;
356
357 /* Decrement the loop counter */
358 i--;
359 } while (i > 0U);
360
361 /* pbuff initialized to input buffer */
362 pbuff = pInlineBuffer;
363
364 /* pS1 initialized to pState */
365 pS1 = pState;
366
367 /* Initializing the loop counter */
368 i = S->N;
369
370 do
371 {
372 /* Writing the re-ordered output back to inplace input buffer */
373 *pbuff++ = *pS1++;
374
375 /* Decrement the loop counter */
376 i--;
377 } while (i > 0U);
378
379
380 /* ---------------------------------------------------------
381 * Step2: Calculate RFFT for N-point input
382 * ---------------------------------------------------------- */
383 /* pInlineBuffer is real input of length N , pState is the complex output of length 2N */
384 arm_rfft_f32 (S->pRfft, pInlineBuffer, pState);
385
386 /*----------------------------------------------------------------------
387 * Step3: Multiply the FFT output with the weights.
388 *----------------------------------------------------------------------*/
389 arm_cmplx_mult_cmplx_f32 (pState, weights, pState, S->N);
390
391 /* ----------- Post-processing ---------- */
392 /* DCT-IV can be obtained from DCT-II by the equation,
393 * Y4(k) = Y2(k) - Y4(k-1) and Y4(-1) = Y4(0)
394 * Hence, Y4(0) = Y2(0)/2 */
395 /* Getting only real part from the output and Converting to DCT-IV */
396
397 /* pbuff initialized to input buffer. */
398 pbuff = pInlineBuffer;
399
400 /* pS1 initialized to pState */
401 pS1 = pState;
402
403 /* Calculating Y4(0) from Y2(0) using Y4(0) = Y2(0)/2 */
404 in = *pS1++ * (float32_t) 0.5;
405 /* input buffer acts as inplace, so output values are stored in the input itself. */
406 *pbuff++ = in;
407
408 /* pState pointer is incremented twice as the real values are located alternatively in the array */
409 pS1++;
410
411 /* Initializing the loop counter */
412 i = (S->N - 1U);
413
414 do
415 {
416 /* Calculating Y4(1) to Y4(N-1) from Y2 using equation Y4(k) = Y2(k) - Y4(k-1) */
417 /* pState pointer (pS1) is incremented twice as the real values are located alternatively in the array */
418 in = *pS1++ - in;
419 *pbuff++ = in;
420
421 /* points to the next real value */
422 pS1++;
423
424 /* Decrement loop counter */
425 i--;
426 } while (i > 0U);
427
428 /*------------ Normalizing the output by multiplying with the normalizing factor ----------*/
429
430 /* Initializing loop counter */
431 i = S->N;
432
433 /* pbuff initialized to the pInlineBuffer (now contains the output values) */
434 pbuff = pInlineBuffer;
435
436 do
437 {
438 /* Multiplying pInlineBuffer with the normalizing factor sqrt(2/N) */
439 in = *pbuff;
440 *pbuff++ = in * S->normalize;
441
442 /* Decrement loop counter */
443 i--;
444 } while (i > 0U);
445
446 #endif /* #if defined (ARM_MATH_LOOPUNROLL) */
447
448 }
449
450 /**
451 @} end of DCT4F32 group
452 */
453
