1 /* ----------------------------------------------------------------------
2  * Project:      CMSIS DSP Library
3  * Title:        arm_dct4_q15.c
4  * Description:  Processing function of DCT4 & IDCT4 Q15
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   @addtogroup DCT4_IDCT4
33   @{
34  */
35 
36 /**
37   @brief         Processing function for the Q15 DCT4/IDCT4.
38   @param[in]     S             points to an instance of the Q15 DCT4 structure.
39   @param[in]     pState        points to state buffer.
40   @param[in,out] pInlineBuffer points to the in-place input and output buffer.
41   @return        none
42 
43   @par           Input an output formats
44                    Internally inputs are downscaled in the RFFT process function to avoid overflows.
45                    Number of bits downscaled, depends on the size of the transform. The input and output
46                    formats for different DCT sizes and number of bits to upscale are mentioned in the table below:
47 
48                    \image html dct4FormatsQ15Table.gif
49  */
50 
arm_dct4_q15(const arm_dct4_instance_q15 * S,q15_t * pState,q15_t * pInlineBuffer)51 void arm_dct4_q15(
52   const arm_dct4_instance_q15 * S,
53         q15_t * pState,
54         q15_t * pInlineBuffer)
55 {
56   const q15_t *weights = S->pTwiddle;                  /* Pointer to the Weights table */
57   const q15_t *cosFact = S->pCosFactor;                /* Pointer to the cos factors table */
58         q15_t *pS1, *pS2, *pbuff;                      /* Temporary pointers for input buffer and pState buffer */
59         q15_t in;                                      /* Temporary variable */
60         uint32_t i;                                    /* Loop counter */
61 
62 
63   /* DCT4 computation involves DCT2 (which is calculated using RFFT)
64    * along with some pre-processing and post-processing.
65    * Computational procedure is explained as follows:
66    * (a) Pre-processing involves multiplying input with cos factor,
67    *     r(n) = 2 * u(n) * cos(pi*(2*n+1)/(4*n))
68    *              where,
69    *                 r(n) -- output of preprocessing
70    *                 u(n) -- input to preprocessing(actual Source buffer)
71    * (b) Calculation of DCT2 using FFT is divided into three steps:
72    *                  Step1: Re-ordering of even and odd elements of input.
73    *                  Step2: Calculating FFT of the re-ordered input.
74    *                  Step3: Taking the real part of the product of FFT output and weights.
75    * (c) Post-processing - DCT4 can be obtained from DCT2 output using the following equation:
76    *                   Y4(k) = Y2(k) - Y4(k-1) and Y4(-1) = Y4(0)
77    *                        where,
78    *                           Y4 -- DCT4 output,   Y2 -- DCT2 output
79    * (d) Multiplying the output with the normalizing factor sqrt(2/N).
80    */
81 
82   /*-------- Pre-processing ------------*/
83   /* Multiplying input with cos factor i.e. r(n) = 2 * x(n) * cos(pi*(2*n+1)/(4*n)) */
84   arm_mult_q15 (pInlineBuffer, cosFact, pInlineBuffer, S->N);
85   arm_shift_q15 (pInlineBuffer, 1, pInlineBuffer, S->N);
86 
87   /* ----------------------------------------------------------------
88    * Step1: Re-ordering of even and odd elements as
89    *             pState[i] =  pInlineBuffer[2*i] and
90    *             pState[N-i-1] = pInlineBuffer[2*i+1] where i = 0 to N/2
91    ---------------------------------------------------------------------*/
92 
93   /* pS1 initialized to pState */
94   pS1 = pState;
95 
96   /* pS2 initialized to pState+N-1, so that it points to the end of the state buffer */
97   pS2 = pState + (S->N - 1U);
98 
99   /* pbuff initialized to input buffer */
100   pbuff = pInlineBuffer;
101 
102 
103 #if defined (ARM_MATH_LOOPUNROLL)
104 
105   /* Initializing the loop counter to N/2 >> 2 for loop unrolling by 4 */
106   i = S->Nby2 >> 2U;
107 
108   /* First part of the processing with loop unrolling.  Compute 4 outputs at a time.
109    ** a second loop below computes the remaining 1 to 3 samples. */
110   do
111   {
112     /* Re-ordering of even and odd elements */
113     /* pState[i] =  pInlineBuffer[2*i] */
114     *pS1++ = *pbuff++;
115     /* pState[N-i-1] = pInlineBuffer[2*i+1] */
116     *pS2-- = *pbuff++;
117 
118     *pS1++ = *pbuff++;
119     *pS2-- = *pbuff++;
120 
121     *pS1++ = *pbuff++;
122     *pS2-- = *pbuff++;
123 
124     *pS1++ = *pbuff++;
125     *pS2-- = *pbuff++;
126 
127     /* Decrement loop counter */
128     i--;
129   } while (i > 0U);
130 
131   /* pbuff initialized to input buffer */
132   pbuff = pInlineBuffer;
133 
134   /* pS1 initialized to pState */
135   pS1 = pState;
136 
137   /* Initializing the loop counter to N/4 instead of N for loop unrolling */
138   i = S->N >> 2U;
139 
140   /* Processing with loop unrolling 4 times as N is always multiple of 4.
141    * Compute 4 outputs at a time */
142   do
143   {
144     /* Writing the re-ordered output back to inplace input buffer */
145     *pbuff++ = *pS1++;
146     *pbuff++ = *pS1++;
147     *pbuff++ = *pS1++;
148     *pbuff++ = *pS1++;
149 
150     /* Decrement the loop counter */
151     i--;
152   } while (i > 0U);
153 
154 
155   /* ---------------------------------------------------------
156    *     Step2: Calculate RFFT for N-point input
157    * ---------------------------------------------------------- */
158   /* pInlineBuffer is real input of length N , pState is the complex output of length 2N */
159   arm_rfft_q15 (S->pRfft, pInlineBuffer, pState);
160 
161   /*----------------------------------------------------------------------
162    *  Step3: Multiply the FFT output with the weights.
163    *----------------------------------------------------------------------*/
164   arm_cmplx_mult_cmplx_q15 (pState, weights, pState, S->N);
165 
166   /* The output of complex multiplication is in 3.13 format.
167    * Hence changing the format of N (i.e. 2*N elements) complex numbers to 1.15 format by shifting left by 2 bits. */
168   arm_shift_q15 (pState, 2, pState, S->N * 2);
169 
170   /* ----------- Post-processing ---------- */
171   /* DCT-IV can be obtained from DCT-II by the equation,
172    *       Y4(k) = Y2(k) - Y4(k-1) and Y4(-1) = Y4(0)
173    *       Hence, Y4(0) = Y2(0)/2  */
174   /* Getting only real part from the output and Converting to DCT-IV */
175 
176   /* Initializing the loop counter to N >> 2 for loop unrolling by 4 */
177   i = (S->N - 1U) >> 2U;
178 
179   /* pbuff initialized to input buffer. */
180   pbuff = pInlineBuffer;
181 
182   /* pS1 initialized to pState */
183   pS1 = pState;
184 
185   /* Calculating Y4(0) from Y2(0) using Y4(0) = Y2(0)/2 */
186   in = *pS1++ >> 1U;
187   /* input buffer acts as inplace, so output values are stored in the input itself. */
188   *pbuff++ = in;
189 
190   /* pState pointer is incremented twice as the real values are located alternatively in the array */
191   pS1++;
192 
193   /* First part of the processing with loop unrolling.  Compute 4 outputs at a time.
194    ** a second loop below computes the remaining 1 to 3 samples. */
195   do
196   {
197     /* Calculating Y4(1) to Y4(N-1) from Y2 using equation Y4(k) = Y2(k) - Y4(k-1) */
198     /* pState pointer (pS1) is incremented twice as the real values are located alternatively in the array */
199     in = *pS1++ - in;
200     *pbuff++ = in;
201     /* points to the next real value */
202     pS1++;
203 
204     in = *pS1++ - in;
205     *pbuff++ = in;
206     pS1++;
207 
208     in = *pS1++ - in;
209     *pbuff++ = in;
210     pS1++;
211 
212     in = *pS1++ - in;
213     *pbuff++ = in;
214     pS1++;
215 
216     /* Decrement the loop counter */
217     i--;
218   } while (i > 0U);
219 
220   /* If the blockSize is not a multiple of 4, compute any remaining output samples here.
221    ** No loop unrolling is used. */
222   i = (S->N - 1U) % 0x4U;
223 
224   while (i > 0U)
225   {
226     /* Calculating Y4(1) to Y4(N-1) from Y2 using equation Y4(k) = Y2(k) - Y4(k-1) */
227     /* pState pointer (pS1) is incremented twice as the real values are located alternatively in the array */
228     in = *pS1++ - in;
229     *pbuff++ = in;
230 
231     /* points to the next real value */
232     pS1++;
233 
234     /* Decrement loop counter */
235     i--;
236   }
237 
238 
239   /*------------ Normalizing the output by multiplying with the normalizing factor ----------*/
240 
241   /* Initializing the loop counter to N/4 instead of N for loop unrolling */
242   i = S->N >> 2U;
243 
244   /* pbuff initialized to the pInlineBuffer(now contains the output values) */
245   pbuff = pInlineBuffer;
246 
247   /* Processing with loop unrolling 4 times as N is always multiple of 4.  Compute 4 outputs at a time */
248   do
249   {
250     /* Multiplying pInlineBuffer with the normalizing factor sqrt(2/N) */
251     in = *pbuff;
252     *pbuff++ = ((q15_t) (((q31_t) in * S->normalize) >> 15));
253 
254     in = *pbuff;
255     *pbuff++ = ((q15_t) (((q31_t) in * S->normalize) >> 15));
256 
257     in = *pbuff;
258     *pbuff++ = ((q15_t) (((q31_t) in * S->normalize) >> 15));
259 
260     in = *pbuff;
261     *pbuff++ = ((q15_t) (((q31_t) in * S->normalize) >> 15));
262 
263     /* Decrement loop counter */
264     i--;
265   } while (i > 0U);
266 
267 
268 #else
269 
270   /* Initializing the loop counter to N/2 */
271   i = S->Nby2;
272 
273   do
274   {
275     /* Re-ordering of even and odd elements */
276     /* pState[i] =  pInlineBuffer[2*i] */
277     *pS1++ = *pbuff++;
278     /* pState[N-i-1] = pInlineBuffer[2*i+1] */
279     *pS2-- = *pbuff++;
280 
281     /* Decrement the loop counter */
282     i--;
283   } while (i > 0U);
284 
285   /* pbuff initialized to input buffer */
286   pbuff = pInlineBuffer;
287 
288   /* pS1 initialized to pState */
289   pS1 = pState;
290 
291   /* Initializing the loop counter */
292   i = S->N;
293 
294   do
295   {
296     /* Writing the re-ordered output back to inplace input buffer */
297     *pbuff++ = *pS1++;
298 
299     /* Decrement the loop counter */
300     i--;
301   } while (i > 0U);
302 
303 
304   /* ---------------------------------------------------------
305    *     Step2: Calculate RFFT for N-point input
306    * ---------------------------------------------------------- */
307   /* pInlineBuffer is real input of length N , pState is the complex output of length 2N */
308   arm_rfft_q15 (S->pRfft, pInlineBuffer, pState);
309 
310   /*----------------------------------------------------------------------
311    *  Step3: Multiply the FFT output with the weights.
312    *----------------------------------------------------------------------*/
313   arm_cmplx_mult_cmplx_q15 (pState, weights, pState, S->N);
314 
315   /* The output of complex multiplication is in 3.13 format.
316    * Hence changing the format of N (i.e. 2*N elements) complex numbers to 1.15 format by shifting left by 2 bits. */
317   arm_shift_q15 (pState, 2, pState, S->N * 2);
318 
319   /* ----------- Post-processing ---------- */
320   /* DCT-IV can be obtained from DCT-II by the equation,
321    *       Y4(k) = Y2(k) - Y4(k-1) and Y4(-1) = Y4(0)
322    *       Hence, Y4(0) = Y2(0)/2  */
323   /* Getting only real part from the output and Converting to DCT-IV */
324 
325   /* pbuff initialized to input buffer. */
326   pbuff = pInlineBuffer;
327 
328   /* pS1 initialized to pState */
329   pS1 = pState;
330 
331   /* Calculating Y4(0) from Y2(0) using Y4(0) = Y2(0)/2 */
332   in = *pS1++ >> 1U;
333   /* input buffer acts as inplace, so output values are stored in the input itself. */
334   *pbuff++ = in;
335 
336   /* pState pointer is incremented twice as the real values are located alternatively in the array */
337   pS1++;
338 
339   /* Initializing the loop counter */
340   i = (S->N - 1U);
341 
342   do
343   {
344     /* Calculating Y4(1) to Y4(N-1) from Y2 using equation Y4(k) = Y2(k) - Y4(k-1) */
345     /* pState pointer (pS1) is incremented twice as the real values are located alternatively in the array */
346     in = *pS1++ - in;
347     *pbuff++ = in;
348 
349     /* points to the next real value */
350     pS1++;
351 
352     /* Decrement loop counter */
353     i--;
354   } while (i > 0U);
355 
356   /*------------ Normalizing the output by multiplying with the normalizing factor ----------*/
357 
358   /* Initializing loop counter */
359   i = S->N;
360 
361   /* pbuff initialized to the pInlineBuffer (now contains the output values) */
362   pbuff = pInlineBuffer;
363 
364   do
365   {
366     /* Multiplying pInlineBuffer with the normalizing factor sqrt(2/N) */
367     in = *pbuff;
368     *pbuff++ = ((q15_t) (((q31_t) in * S->normalize) >> 15));
369 
370     /* Decrement loop counter */
371     i--;
372 
373   } while (i > 0U);
374 
375 #endif /* #if defined (ARM_MATH_LOOPUNROLL) */
376 
377 }
378 
379 /**
380   @} end of DCT4_IDCT4 group
381  */
382