1 /*
2  * Copyright (C) 2015 The Android Open Source Project
3  *
4  * Licensed under the Apache License, Version 2.0 (the "License");
5  * you may not use this file except in compliance with the License.
6  * You may obtain a copy of the License at
7  *
8  *      http://www.apache.org/licenses/LICENSE-2.0
9  *
10  * Unless required by applicable law or agreed to in writing, software
11  * distributed under the License is distributed on an "AS IS" BASIS,
12  * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
13  * See the License for the specific language governing permissions and
14  * limitations under the License.
15  */
16 
17 package com.android.cts.verifier.audio.wavelib;
18 
19 public class DspFftServer {
20     private int mN = 0;
21     private int mOrder = 0;
22 
23     DspBufferDouble mCos;
24     DspBufferDouble mSin;
25     public boolean isInitialized = false;
26 
DspFftServer(int size)27     public DspFftServer(int size) {
28         init(size);
29     }
30 
init(int size)31     public boolean init(int size) {
32         boolean status = false;
33         mN=size;
34 
35         mOrder = (int) (Math.log(mN) / Math.log(2));
36         if (mN == (1 << mOrder)) {
37             mCos = new DspBufferDouble(mN / 2);
38             mSin = new DspBufferDouble(mN / 2);
39             for (int i = 0; i < mN / 2; i++) {
40                 mCos.mData[i] = Math.cos(-2 * Math.PI * i / mN);
41                 mSin.mData[i] = Math.sin(-2 * Math.PI * i / mN);
42             }
43             status = true;
44         } else {
45             mN = 0;
46             throw new RuntimeException("FFT must be power of 2");
47         }
48         isInitialized = status;
49         return status;
50     }
51 
fft(DspBufferComplex r, int sign)52     public void fft(DspBufferComplex r, int sign) {
53         int ii, jj, kk, n1, n2, aa;
54         double cc, ss, t1, t2;
55 
56         // Bit-reverse
57         jj = 0;
58         n2 = mN / 2;
59         for (ii = 1; ii < mN - 1; ii++) {
60             n1 = n2;
61             while (jj >= n1) {
62                 jj = jj - n1;
63                 n1 = n1 / 2;
64             }
65             jj = jj + n1;
66 
67             if (ii < jj) {
68                 t1 =  r.mReal[ii];
69                 r.mReal[ii] = r.mReal[jj];
70                 r.mReal[jj] = t1;
71                 t1 = r.mImag[ii];
72                 r.mImag[ii] = r.mImag[jj];
73                 r.mImag[jj] = t1;
74             }
75         }
76 
77         // FFT
78         n1 = 0;
79         n2 = 1;
80         for (ii = 0; ii < mOrder; ii++) {
81             n1 = n2;
82             n2 = n2 + n2;
83             aa = 0;
84 
85             for (jj = 0; jj < n1; jj++) {
86                 cc = mCos.mData[aa];
87                 ss = sign * mSin.mData[aa];
88                 aa += 1 << (mOrder - ii - 1);
89                 for (kk = jj; kk < mN; kk = kk + n2) {
90                     t1 = cc * r.mReal[kk + n1] - ss * r.mImag[kk + n1];
91                     t2 = ss * r.mReal[kk + n1] + cc * r.mImag[kk + n1];
92                     r.mReal[kk + n1] = r.mReal[kk] - t1;
93                     r.mImag[kk + n1] = r.mImag[kk] - t2;
94                     r.mReal[kk] = r.mReal[kk] + t1;
95                     r.mImag[kk] = r.mImag[kk] + t2;
96                 }
97             }
98         }
99     }
100 }
101