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 @return none
130 */
131
arm_dct4_f32(const arm_dct4_instance_f32 * S,float32_t * pState,float32_t * pInlineBuffer)132 void arm_dct4_f32(
133 const arm_dct4_instance_f32 * S,
134 float32_t * pState,
135 float32_t * pInlineBuffer)
136 {
137 const float32_t *weights = S->pTwiddle; /* Pointer to the Weights table */
138 const float32_t *cosFact = S->pCosFactor; /* Pointer to the cos factors table */
139 float32_t *pS1, *pS2, *pbuff; /* Temporary pointers for input buffer and pState buffer */
140 float32_t in; /* Temporary variable */
141 uint32_t i; /* Loop counter */
142
143
144 /* DCT4 computation involves DCT2 (which is calculated using RFFT)
145 * along with some pre-processing and post-processing.
146 * Computational procedure is explained as follows:
147 * (a) Pre-processing involves multiplying input with cos factor,
148 * r(n) = 2 * u(n) * cos(pi*(2*n+1)/(4*n))
149 * where,
150 * r(n) -- output of preprocessing
151 * u(n) -- input to preprocessing(actual Source buffer)
152 * (b) Calculation of DCT2 using FFT is divided into three steps:
153 * Step1: Re-ordering of even and odd elements of input.
154 * Step2: Calculating FFT of the re-ordered input.
155 * Step3: Taking the real part of the product of FFT output and weights.
156 * (c) Post-processing - DCT4 can be obtained from DCT2 output using the following equation:
157 * Y4(k) = Y2(k) - Y4(k-1) and Y4(-1) = Y4(0)
158 * where,
159 * Y4 -- DCT4 output, Y2 -- DCT2 output
160 * (d) Multiplying the output with the normalizing factor sqrt(2/N).
161 */
162
163 /*-------- Pre-processing ------------*/
164 /* Multiplying input with cos factor i.e. r(n) = 2 * x(n) * cos(pi*(2*n+1)/(4*n)) */
165 arm_scale_f32(pInlineBuffer, 2.0f, pInlineBuffer, S->N);
166 arm_mult_f32(pInlineBuffer, cosFact, pInlineBuffer, S->N);
167
168 /* ----------------------------------------------------------------
169 * Step1: Re-ordering of even and odd elements as
170 * pState[i] = pInlineBuffer[2*i] and
171 * pState[N-i-1] = pInlineBuffer[2*i+1] where i = 0 to N/2
172 ---------------------------------------------------------------------*/
173
174 /* pS1 initialized to pState */
175 pS1 = pState;
176
177 /* pS2 initialized to pState+N-1, so that it points to the end of the state buffer */
178 pS2 = pState + (S->N - 1U);
179
180 /* pbuff initialized to input buffer */
181 pbuff = pInlineBuffer;
182
183
184 #if defined (ARM_MATH_LOOPUNROLL)
185
186 /* Initializing the loop counter to N/2 >> 2 for loop unrolling by 4 */
187 i = S->Nby2 >> 2U;
188
189 /* First part of the processing with loop unrolling. Compute 4 outputs at a time.
190 ** a second loop below computes the remaining 1 to 3 samples. */
191 do
192 {
193 /* Re-ordering of even and odd elements */
194 /* pState[i] = pInlineBuffer[2*i] */
195 *pS1++ = *pbuff++;
196 /* pState[N-i-1] = pInlineBuffer[2*i+1] */
197 *pS2-- = *pbuff++;
198
199 *pS1++ = *pbuff++;
200 *pS2-- = *pbuff++;
201
202 *pS1++ = *pbuff++;
203 *pS2-- = *pbuff++;
204
205 *pS1++ = *pbuff++;
206 *pS2-- = *pbuff++;
207
208 /* Decrement loop counter */
209 i--;
210 } while (i > 0U);
211
212 /* pbuff initialized to input buffer */
213 pbuff = pInlineBuffer;
214
215 /* pS1 initialized to pState */
216 pS1 = pState;
217
218 /* Initializing the loop counter to N/4 instead of N for loop unrolling */
219 i = S->N >> 2U;
220
221 /* Processing with loop unrolling 4 times as N is always multiple of 4.
222 * Compute 4 outputs at a time */
223 do
224 {
225 /* Writing the re-ordered output back to inplace input buffer */
226 *pbuff++ = *pS1++;
227 *pbuff++ = *pS1++;
228 *pbuff++ = *pS1++;
229 *pbuff++ = *pS1++;
230
231 /* Decrement the loop counter */
232 i--;
233 } while (i > 0U);
234
235
236 /* ---------------------------------------------------------
237 * Step2: Calculate RFFT for N-point input
238 * ---------------------------------------------------------- */
239 /* pInlineBuffer is real input of length N , pState is the complex output of length 2N */
240 arm_rfft_f32 (S->pRfft, pInlineBuffer, pState);
241
242 /*----------------------------------------------------------------------
243 * Step3: Multiply the FFT output with the weights.
244 *----------------------------------------------------------------------*/
245 arm_cmplx_mult_cmplx_f32 (pState, weights, pState, S->N);
246
247 /* ----------- Post-processing ---------- */
248 /* DCT-IV can be obtained from DCT-II by the equation,
249 * Y4(k) = Y2(k) - Y4(k-1) and Y4(-1) = Y4(0)
250 * Hence, Y4(0) = Y2(0)/2 */
251 /* Getting only real part from the output and Converting to DCT-IV */
252
253 /* Initializing the loop counter to N >> 2 for loop unrolling by 4 */
254 i = (S->N - 1U) >> 2U;
255
256 /* pbuff initialized to input buffer. */
257 pbuff = pInlineBuffer;
258
259 /* pS1 initialized to pState */
260 pS1 = pState;
261
262 /* Calculating Y4(0) from Y2(0) using Y4(0) = Y2(0)/2 */
263 in = *pS1++ * (float32_t) 0.5;
264 /* input buffer acts as inplace, so output values are stored in the input itself. */
265 *pbuff++ = in;
266
267 /* pState pointer is incremented twice as the real values are located alternatively in the array */
268 pS1++;
269
270 /* First part of the processing with loop unrolling. Compute 4 outputs at a time.
271 ** a second loop below computes the remaining 1 to 3 samples. */
272 do
273 {
274 /* Calculating Y4(1) to Y4(N-1) from Y2 using equation Y4(k) = Y2(k) - Y4(k-1) */
275 /* pState pointer (pS1) is incremented twice as the real values are located alternatively in the array */
276 in = *pS1++ - in;
277 *pbuff++ = in;
278 /* points to the next real value */
279 pS1++;
280
281 in = *pS1++ - in;
282 *pbuff++ = in;
283 pS1++;
284
285 in = *pS1++ - in;
286 *pbuff++ = in;
287 pS1++;
288
289 in = *pS1++ - in;
290 *pbuff++ = in;
291 pS1++;
292
293 /* Decrement the loop counter */
294 i--;
295 } while (i > 0U);
296
297 /* If the blockSize is not a multiple of 4, compute any remaining output samples here.
298 ** No loop unrolling is used. */
299 i = (S->N - 1U) % 0x4U;
300
301 while (i > 0U)
302 {
303 /* Calculating Y4(1) to Y4(N-1) from Y2 using equation Y4(k) = Y2(k) - Y4(k-1) */
304 /* pState pointer (pS1) is incremented twice as the real values are located alternatively in the array */
305 in = *pS1++ - in;
306 *pbuff++ = in;
307
308 /* points to the next real value */
309 pS1++;
310
311 /* Decrement the loop counter */
312 i--;
313 }
314
315
316 /*------------ Normalizing the output by multiplying with the normalizing factor ----------*/
317
318 /* Initializing the loop counter to N/4 instead of N for loop unrolling */
319 i = S->N >> 2U;
320
321 /* pbuff initialized to the pInlineBuffer(now contains the output values) */
322 pbuff = pInlineBuffer;
323
324 /* Processing with loop unrolling 4 times as N is always multiple of 4. Compute 4 outputs at a time */
325 do
326 {
327 /* Multiplying pInlineBuffer with the normalizing factor sqrt(2/N) */
328 in = *pbuff;
329 *pbuff++ = in * S->normalize;
330
331 in = *pbuff;
332 *pbuff++ = in * S->normalize;
333
334 in = *pbuff;
335 *pbuff++ = in * S->normalize;
336
337 in = *pbuff;
338 *pbuff++ = in * S->normalize;
339
340 /* Decrement the loop counter */
341 i--;
342 } while (i > 0U);
343
344
345 #else
346
347 /* Initializing the loop counter to N/2 */
348 i = S->Nby2;
349
350 do
351 {
352 /* Re-ordering of even and odd elements */
353 /* pState[i] = pInlineBuffer[2*i] */
354 *pS1++ = *pbuff++;
355 /* pState[N-i-1] = pInlineBuffer[2*i+1] */
356 *pS2-- = *pbuff++;
357
358 /* Decrement the loop counter */
359 i--;
360 } while (i > 0U);
361
362 /* pbuff initialized to input buffer */
363 pbuff = pInlineBuffer;
364
365 /* pS1 initialized to pState */
366 pS1 = pState;
367
368 /* Initializing the loop counter */
369 i = S->N;
370
371 do
372 {
373 /* Writing the re-ordered output back to inplace input buffer */
374 *pbuff++ = *pS1++;
375
376 /* Decrement the loop counter */
377 i--;
378 } while (i > 0U);
379
380
381 /* ---------------------------------------------------------
382 * Step2: Calculate RFFT for N-point input
383 * ---------------------------------------------------------- */
384 /* pInlineBuffer is real input of length N , pState is the complex output of length 2N */
385 arm_rfft_f32 (S->pRfft, pInlineBuffer, pState);
386
387 /*----------------------------------------------------------------------
388 * Step3: Multiply the FFT output with the weights.
389 *----------------------------------------------------------------------*/
390 arm_cmplx_mult_cmplx_f32 (pState, weights, pState, S->N);
391
392 /* ----------- Post-processing ---------- */
393 /* DCT-IV can be obtained from DCT-II by the equation,
394 * Y4(k) = Y2(k) - Y4(k-1) and Y4(-1) = Y4(0)
395 * Hence, Y4(0) = Y2(0)/2 */
396 /* Getting only real part from the output and Converting to DCT-IV */
397
398 /* pbuff initialized to input buffer. */
399 pbuff = pInlineBuffer;
400
401 /* pS1 initialized to pState */
402 pS1 = pState;
403
404 /* Calculating Y4(0) from Y2(0) using Y4(0) = Y2(0)/2 */
405 in = *pS1++ * (float32_t) 0.5;
406 /* input buffer acts as inplace, so output values are stored in the input itself. */
407 *pbuff++ = in;
408
409 /* pState pointer is incremented twice as the real values are located alternatively in the array */
410 pS1++;
411
412 /* Initializing the loop counter */
413 i = (S->N - 1U);
414
415 do
416 {
417 /* Calculating Y4(1) to Y4(N-1) from Y2 using equation Y4(k) = Y2(k) - Y4(k-1) */
418 /* pState pointer (pS1) is incremented twice as the real values are located alternatively in the array */
419 in = *pS1++ - in;
420 *pbuff++ = in;
421
422 /* points to the next real value */
423 pS1++;
424
425 /* Decrement loop counter */
426 i--;
427 } while (i > 0U);
428
429 /*------------ Normalizing the output by multiplying with the normalizing factor ----------*/
430
431 /* Initializing loop counter */
432 i = S->N;
433
434 /* pbuff initialized to the pInlineBuffer (now contains the output values) */
435 pbuff = pInlineBuffer;
436
437 do
438 {
439 /* Multiplying pInlineBuffer with the normalizing factor sqrt(2/N) */
440 in = *pbuff;
441 *pbuff++ = in * S->normalize;
442
443 /* Decrement loop counter */
444 i--;
445 } while (i > 0U);
446
447 #endif /* #if defined (ARM_MATH_LOOPUNROLL) */
448
449 }
450
451 /**
452 @} end of DCT4F32 group
453 */
454