1 /* Copyright (c) 2008 Xiph.Org Foundation
2 Written by Jean-Marc Valin */
3 /*
4 Redistribution and use in source and binary forms, with or without
5 modification, are permitted provided that the following conditions
6 are met:
7
8 - Redistributions of source code must retain the above copyright
9 notice, this list of conditions and the following disclaimer.
10
11 - Redistributions in binary form must reproduce the above copyright
12 notice, this list of conditions and the following disclaimer in the
13 documentation and/or other materials provided with the distribution.
14
15 THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
16 ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
17 LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
18 A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER
19 OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
20 EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
21 PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
22 PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
23 LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
24 NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
25 SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
26 */
27
28 #ifdef HAVE_CONFIG_H
29 #include "config.h"
30 #endif
31
32 #include <stdio.h>
33
34 #include "stack_alloc.h"
35 #include "kiss_fft.h"
36 #include "mathops.h"
37 #include "modes.h"
38
39 #ifndef M_PI
40 #define M_PI 3.141592653
41 #endif
42
43 int ret = 0;
44
check(kiss_fft_cpx * in,kiss_fft_cpx * out,int nfft,int isinverse)45 void check(kiss_fft_cpx * in,kiss_fft_cpx * out,int nfft,int isinverse)
46 {
47 int bin,k;
48 double errpow=0,sigpow=0, snr;
49
50 for (bin=0;bin<nfft;++bin) {
51 double ansr = 0;
52 double ansi = 0;
53 double difr;
54 double difi;
55
56 for (k=0;k<nfft;++k) {
57 double phase = -2*M_PI*bin*k/nfft;
58 double re = cos(phase);
59 double im = sin(phase);
60 if (isinverse)
61 im = -im;
62
63 if (!isinverse)
64 {
65 re /= nfft;
66 im /= nfft;
67 }
68
69 ansr += in[k].r * re - in[k].i * im;
70 ansi += in[k].r * im + in[k].i * re;
71 }
72 /*printf ("%d %d ", (int)ansr, (int)ansi);*/
73 difr = ansr - out[bin].r;
74 difi = ansi - out[bin].i;
75 errpow += difr*difr + difi*difi;
76 sigpow += ansr*ansr+ansi*ansi;
77 }
78 snr = 10*log10(sigpow/errpow);
79 printf("nfft=%d inverse=%d,snr = %f\n",nfft,isinverse,snr );
80 if (snr<60) {
81 printf( "** poor snr: %f ** \n", snr);
82 ret = 1;
83 }
84 }
85
test1d(int nfft,int isinverse,int arch)86 void test1d(int nfft,int isinverse,int arch)
87 {
88 size_t buflen = sizeof(kiss_fft_cpx)*nfft;
89 kiss_fft_cpx *in;
90 kiss_fft_cpx *out;
91 int k;
92 #ifdef CUSTOM_MODES
93 kiss_fft_state *cfg = opus_fft_alloc(nfft,0,0,arch);
94 #else
95 int id;
96 const kiss_fft_state *cfg;
97 CELTMode *mode = opus_custom_mode_create(48000, 960, NULL);
98 if (nfft == 480) id = 0;
99 else if (nfft == 240) id = 1;
100 else if (nfft == 120) id = 2;
101 else if (nfft == 60) id = 3;
102 else return;
103 cfg = mode->mdct.kfft[id];
104 #endif
105
106 in = (kiss_fft_cpx*)malloc(buflen);
107 out = (kiss_fft_cpx*)malloc(buflen);
108
109 for (k=0;k<nfft;++k) {
110 in[k].r = (rand() % 32767) - 16384;
111 in[k].i = (rand() % 32767) - 16384;
112 }
113
114 for (k=0;k<nfft;++k) {
115 in[k].r *= 32768;
116 in[k].i *= 32768;
117 }
118
119 if (isinverse)
120 {
121 for (k=0;k<nfft;++k) {
122 in[k].r /= nfft;
123 in[k].i /= nfft;
124 }
125 }
126
127 /*for (k=0;k<nfft;++k) printf("%d %d ", in[k].r, in[k].i);printf("\n");*/
128
129 if (isinverse)
130 opus_ifft(cfg,in,out, arch);
131 else
132 opus_fft(cfg,in,out, arch);
133
134 /*for (k=0;k<nfft;++k) printf("%d %d ", out[k].r, out[k].i);printf("\n");*/
135
136 check(in,out,nfft,isinverse);
137
138 free(in);
139 free(out);
140 #ifdef CUSTOM_MODES
141 opus_fft_free(cfg, arch);
142 #endif
143 }
144
main(int argc,char ** argv)145 int main(int argc,char ** argv)
146 {
147 ALLOC_STACK;
148 int arch = opus_select_arch();
149
150 if (argc>1) {
151 int k;
152 for (k=1;k<argc;++k) {
153 test1d(atoi(argv[k]),0,arch);
154 test1d(atoi(argv[k]),1,arch);
155 }
156 }else{
157 test1d(32,0,arch);
158 test1d(32,1,arch);
159 test1d(128,0,arch);
160 test1d(128,1,arch);
161 test1d(256,0,arch);
162 test1d(256,1,arch);
163 #ifndef RADIX_TWO_ONLY
164 test1d(36,0,arch);
165 test1d(36,1,arch);
166 test1d(50,0,arch);
167 test1d(50,1,arch);
168 test1d(60,0,arch);
169 test1d(60,1,arch);
170 test1d(120,0,arch);
171 test1d(120,1,arch);
172 test1d(240,0,arch);
173 test1d(240,1,arch);
174 test1d(480,0,arch);
175 test1d(480,1,arch);
176 #endif
177 }
178 return ret;
179 }
180