Fix for an assertion when running the fixed point tests
[opus.git] / celt / tests / test_unit_dft.c
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
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
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
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 }