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#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 "stack_alloc.h"
42#include "kiss_fft.h"
43#include "kiss_fft.c"
44#include "mathops.c"
45#include "entcode.c"
46
47
48#ifndef M_PI
49#define M_PI 3.141592653
50#endif
51
52int ret = 0;
53
54void check(kiss_fft_cpx  * in,kiss_fft_cpx  * out,int nfft,int isinverse)
55{
56    int bin,k;
57    double errpow=0,sigpow=0, snr;
58
59    for (bin=0;bin<nfft;++bin) {
60        double ansr = 0;
61        double ansi = 0;
62        double difr;
63        double difi;
64
65        for (k=0;k<nfft;++k) {
66            double phase = -2*M_PI*bin*k/nfft;
67            double re = cos(phase);
68            double im = sin(phase);
69            if (isinverse)
70                im = -im;
71
72            if (!isinverse)
73            {
74               re /= nfft;
75               im /= nfft;
76            }
77
78            ansr += in[k].r * re - in[k].i * im;
79            ansi += in[k].r * im + in[k].i * re;
80        }
81        /*printf ("%d %d ", (int)ansr, (int)ansi);*/
82        difr = ansr - out[bin].r;
83        difi = ansi - out[bin].i;
84        errpow += difr*difr + difi*difi;
85        sigpow += ansr*ansr+ansi*ansi;
86    }
87    snr = 10*log10(sigpow/errpow);
88    printf("nfft=%d inverse=%d,snr = %f\n",nfft,isinverse,snr );
89    if (snr<60) {
90       printf( "** poor snr: %f ** \n", snr);
91       ret = 1;
92    }
93}
94
95void test1d(int nfft,int isinverse)
96{
97    size_t buflen = sizeof(kiss_fft_cpx)*nfft;
98
99    kiss_fft_cpx  * in = (kiss_fft_cpx*)malloc(buflen);
100    kiss_fft_cpx  * out= (kiss_fft_cpx*)malloc(buflen);
101    kiss_fft_state *cfg = opus_fft_alloc(nfft,0,0);
102    int k;
103
104    for (k=0;k<nfft;++k) {
105        in[k].r = (rand() % 32767) - 16384;
106        in[k].i = (rand() % 32767) - 16384;
107    }
108
109    for (k=0;k<nfft;++k) {
110       in[k].r *= 32768;
111       in[k].i *= 32768;
112    }
113
114    if (isinverse)
115    {
116       for (k=0;k<nfft;++k) {
117          in[k].r /= nfft;
118          in[k].i /= nfft;
119       }
120    }
121
122    /*for (k=0;k<nfft;++k) printf("%d %d ", in[k].r, in[k].i);printf("\n");*/
123
124    if (isinverse)
125       opus_ifft(cfg,in,out);
126    else
127       opus_fft(cfg,in,out);
128
129    /*for (k=0;k<nfft;++k) printf("%d %d ", out[k].r, out[k].i);printf("\n");*/
130
131    check(in,out,nfft,isinverse);
132
133    free(in);
134    free(out);
135    free(cfg);
136}
137
138int main(int argc,char ** argv)
139{
140    ALLOC_STACK;
141    if (argc>1) {
142        int k;
143        for (k=1;k<argc;++k) {
144            test1d(atoi(argv[k]),0);
145            test1d(atoi(argv[k]),1);
146        }
147    }else{
148        test1d(32,0);
149        test1d(32,1);
150        test1d(128,0);
151        test1d(128,1);
152        test1d(256,0);
153        test1d(256,1);
154#ifndef RADIX_TWO_ONLY
155        test1d(36,0);
156        test1d(36,1);
157        test1d(50,0);
158        test1d(50,1);
159        test1d(120,0);
160        test1d(120,1);
161#endif
162    }
163    return ret;
164}
165