1/*
2 *  Copyright (c) 2011 The WebRTC project authors. All Rights Reserved.
3 *
4 *  Use of this source code is governed by a BSD-style license
5 *  that can be found in the LICENSE file in the root of the source
6 *  tree. An additional intellectual property rights grant can be found
7 *  in the file PATENTS.  All contributing project authors may
8 *  be found in the AUTHORS file in the root of the source tree.
9 */
10
11#include "settings.h"
12#include "fft.h"
13#include "codec.h"
14#include "os_specific_inline.h"
15#include <math.h>
16
17static double costab1[FRAMESAMPLES_HALF];
18static double sintab1[FRAMESAMPLES_HALF];
19static double costab2[FRAMESAMPLES_QUARTER];
20static double sintab2[FRAMESAMPLES_QUARTER];
21
22void WebRtcIsac_InitTransform()
23{
24  int k;
25  double fact, phase;
26
27  fact = PI / (FRAMESAMPLES_HALF);
28  phase = 0.0;
29  for (k = 0; k < FRAMESAMPLES_HALF; k++) {
30    costab1[k] = cos(phase);
31    sintab1[k] = sin(phase);
32    phase += fact;
33  }
34
35  fact = PI * ((double) (FRAMESAMPLES_HALF - 1)) / ((double) FRAMESAMPLES_HALF);
36  phase = 0.5 * fact;
37  for (k = 0; k < FRAMESAMPLES_QUARTER; k++) {
38    costab2[k] = cos(phase);
39    sintab2[k] = sin(phase);
40    phase += fact;
41  }
42}
43
44
45void WebRtcIsac_Time2Spec(double *inre1,
46                         double *inre2,
47                         int16_t *outreQ7,
48                         int16_t *outimQ7,
49                         FFTstr *fftstr_obj)
50{
51
52  int k;
53  int dims[1];
54  double tmp1r, tmp1i, xr, xi, yr, yi, fact;
55  double tmpre[FRAMESAMPLES_HALF], tmpim[FRAMESAMPLES_HALF];
56
57
58  dims[0] = FRAMESAMPLES_HALF;
59
60
61  /* Multiply with complex exponentials and combine into one complex vector */
62  fact = 0.5 / sqrt(FRAMESAMPLES_HALF);
63  for (k = 0; k < FRAMESAMPLES_HALF; k++) {
64    tmp1r = costab1[k];
65    tmp1i = sintab1[k];
66    tmpre[k] = (inre1[k] * tmp1r + inre2[k] * tmp1i) * fact;
67    tmpim[k] = (inre2[k] * tmp1r - inre1[k] * tmp1i) * fact;
68  }
69
70
71  /* Get DFT */
72  WebRtcIsac_Fftns(1, dims, tmpre, tmpim, -1, 1.0, fftstr_obj);
73
74  /* Use symmetry to separate into two complex vectors and center frames in time around zero */
75  for (k = 0; k < FRAMESAMPLES_QUARTER; k++) {
76    xr = tmpre[k] + tmpre[FRAMESAMPLES_HALF - 1 - k];
77    yi = -tmpre[k] + tmpre[FRAMESAMPLES_HALF - 1 - k];
78    xi = tmpim[k] - tmpim[FRAMESAMPLES_HALF - 1 - k];
79    yr = tmpim[k] + tmpim[FRAMESAMPLES_HALF - 1 - k];
80
81    tmp1r = costab2[k];
82    tmp1i = sintab2[k];
83    outreQ7[k] = (int16_t)WebRtcIsac_lrint((xr * tmp1r - xi * tmp1i) * 128.0);
84    outimQ7[k] = (int16_t)WebRtcIsac_lrint((xr * tmp1i + xi * tmp1r) * 128.0);
85    outreQ7[FRAMESAMPLES_HALF - 1 - k] = (int16_t)WebRtcIsac_lrint((-yr * tmp1i - yi * tmp1r) * 128.0);
86    outimQ7[FRAMESAMPLES_HALF - 1 - k] = (int16_t)WebRtcIsac_lrint((-yr * tmp1r + yi * tmp1i) * 128.0);
87  }
88}
89
90
91void WebRtcIsac_Spec2time(double *inre, double *inim, double *outre1, double *outre2, FFTstr *fftstr_obj)
92{
93
94  int k;
95  double tmp1r, tmp1i, xr, xi, yr, yi, fact;
96
97  int dims;
98
99  dims = FRAMESAMPLES_HALF;
100
101  for (k = 0; k < FRAMESAMPLES_QUARTER; k++) {
102    /* Move zero in time to beginning of frames */
103    tmp1r = costab2[k];
104    tmp1i = sintab2[k];
105    xr = inre[k] * tmp1r + inim[k] * tmp1i;
106    xi = inim[k] * tmp1r - inre[k] * tmp1i;
107    yr = -inim[FRAMESAMPLES_HALF - 1 - k] * tmp1r - inre[FRAMESAMPLES_HALF - 1 - k] * tmp1i;
108    yi = -inre[FRAMESAMPLES_HALF - 1 - k] * tmp1r + inim[FRAMESAMPLES_HALF - 1 - k] * tmp1i;
109
110    /* Combine into one vector,  z = x + j * y */
111    outre1[k] = xr - yi;
112    outre1[FRAMESAMPLES_HALF - 1 - k] = xr + yi;
113    outre2[k] = xi + yr;
114    outre2[FRAMESAMPLES_HALF - 1 - k] = -xi + yr;
115  }
116
117
118  /* Get IDFT */
119  WebRtcIsac_Fftns(1, &dims, outre1, outre2, 1, FRAMESAMPLES_HALF, fftstr_obj);
120
121
122  /* Demodulate and separate */
123  fact = sqrt(FRAMESAMPLES_HALF);
124  for (k = 0; k < FRAMESAMPLES_HALF; k++) {
125    tmp1r = costab1[k];
126    tmp1i = sintab1[k];
127    xr = (outre1[k] * tmp1r - outre2[k] * tmp1i) * fact;
128    outre2[k] = (outre2[k] * tmp1r + outre1[k] * tmp1i) * fact;
129    outre1[k] = xr;
130  }
131}
132