1/* Copyright (c) 2008-2011 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#define SKIP_CONFIG_H
33
34#ifndef CUSTOM_MODES
35#define CUSTOM_MODES
36#endif
37
38#include <stdio.h>
39
40#define CELT_C
41#include "mdct.h"
42#include "stack_alloc.h"
43
44#include "kiss_fft.c"
45#include "mdct.c"
46#include "mathops.c"
47#include "entcode.c"
48
49#ifndef M_PI
50#define M_PI 3.141592653
51#endif
52
53int ret = 0;
54void check(kiss_fft_scalar  * in,kiss_fft_scalar  * out,int nfft,int isinverse)
55{
56    int bin,k;
57    double errpow=0,sigpow=0;
58    double snr;
59    for (bin=0;bin<nfft/2;++bin) {
60        double ansr = 0;
61        double difr;
62
63        for (k=0;k<nfft;++k) {
64           double phase = 2*M_PI*(k+.5+.25*nfft)*(bin+.5)/nfft;
65           double re = cos(phase);
66
67           re /= nfft/4;
68
69           ansr += in[k] * re;
70        }
71        /*printf ("%f %f\n", ansr, out[bin]);*/
72        difr = ansr - out[bin];
73        errpow += difr*difr;
74        sigpow += ansr*ansr;
75    }
76    snr = 10*log10(sigpow/errpow);
77    printf("nfft=%d inverse=%d,snr = %f\n",nfft,isinverse,snr );
78    if (snr<60) {
79       printf( "** poor snr: %f **\n", snr);
80       ret = 1;
81    }
82}
83
84void check_inv(kiss_fft_scalar  * in,kiss_fft_scalar  * out,int nfft,int isinverse)
85{
86   int bin,k;
87   double errpow=0,sigpow=0;
88   double snr;
89   for (bin=0;bin<nfft;++bin) {
90      double ansr = 0;
91      double difr;
92
93      for (k=0;k<nfft/2;++k) {
94         double phase = 2*M_PI*(bin+.5+.25*nfft)*(k+.5)/nfft;
95         double re = cos(phase);
96
97         /*re *= 2;*/
98
99         ansr += in[k] * re;
100      }
101      /*printf ("%f %f\n", ansr, out[bin]);*/
102      difr = ansr - out[bin];
103      errpow += difr*difr;
104      sigpow += ansr*ansr;
105   }
106   snr = 10*log10(sigpow/errpow);
107   printf("nfft=%d inverse=%d,snr = %f\n",nfft,isinverse,snr );
108   if (snr<60) {
109      printf( "** poor snr: %f **\n", snr);
110      ret = 1;
111   }
112}
113
114
115void test1d(int nfft,int isinverse)
116{
117    mdct_lookup cfg;
118    size_t buflen = sizeof(kiss_fft_scalar)*nfft;
119
120    kiss_fft_scalar  * in = (kiss_fft_scalar*)malloc(buflen);
121    kiss_fft_scalar  * in_copy = (kiss_fft_scalar*)malloc(buflen);
122    kiss_fft_scalar  * out= (kiss_fft_scalar*)malloc(buflen);
123    opus_val16  * window= (opus_val16*)malloc(sizeof(opus_val16)*nfft/2);
124    int k;
125
126    clt_mdct_init(&cfg, nfft, 0);
127    for (k=0;k<nfft;++k) {
128        in[k] = (rand() % 32768) - 16384;
129    }
130
131    for (k=0;k<nfft/2;++k) {
132       window[k] = Q15ONE;
133    }
134    for (k=0;k<nfft;++k) {
135       in[k] *= 32768;
136    }
137
138    if (isinverse)
139    {
140       for (k=0;k<nfft;++k) {
141          in[k] /= nfft;
142       }
143    }
144
145    for (k=0;k<nfft;++k)
146       in_copy[k] = in[k];
147    /*for (k=0;k<nfft;++k) printf("%d %d ", in[k].r, in[k].i);printf("\n");*/
148
149    if (isinverse)
150    {
151       for (k=0;k<nfft;++k)
152          out[k] = 0;
153       clt_mdct_backward(&cfg,in,out, window, nfft/2, 0, 1);
154       /* apply TDAC because clt_mdct_backward() no longer does that */
155       for (k=0;k<nfft/4;++k)
156          out[nfft-k-1] = out[nfft/2+k];
157       check_inv(in,out,nfft,isinverse);
158    } else {
159       clt_mdct_forward(&cfg,in,out,window, nfft/2, 0, 1);
160       check(in_copy,out,nfft,isinverse);
161    }
162    /*for (k=0;k<nfft;++k) printf("%d %d ", out[k].r, out[k].i);printf("\n");*/
163
164
165    free(in);
166    free(out);
167    clt_mdct_clear(&cfg);
168}
169
170int main(int argc,char ** argv)
171{
172    ALLOC_STACK;
173    if (argc>1) {
174        int k;
175        for (k=1;k<argc;++k) {
176            test1d(atoi(argv[k]),0);
177            test1d(atoi(argv[k]),1);
178        }
179    }else{
180        test1d(32,0);
181        test1d(32,1);
182        test1d(256,0);
183        test1d(256,1);
184        test1d(512,0);
185        test1d(512,1);
186        test1d(1024,0);
187        test1d(1024,1);
188        test1d(2048,0);
189        test1d(2048,1);
190#ifndef RADIX_TWO_ONLY
191        test1d(36,0);
192        test1d(36,1);
193        test1d(40,0);
194        test1d(40,1);
195        test1d(60,0);
196        test1d(60,1);
197        test1d(120,0);
198        test1d(120,1);
199        test1d(240,0);
200        test1d(240,1);
201        test1d(480,0);
202        test1d(480,1);
203        test1d(960,0);
204        test1d(960,1);
205        test1d(1920,0);
206        test1d(1920,1);
207#endif
208    }
209    return ret;
210}
211