c902769f817501a45b2bd9a90605d6b26ecf4958
[Faustine.git] / interpretor / libsndfile-1.0.25 / tests / dft_cmp.c
1 /*
2 ** Copyright (C) 2002-2011 Erik de Castro Lopo <erikd@mega-nerd.com>
3 **
4 ** This program is free software; you can redistribute it and/or modify
5 ** it under the terms of the GNU General Public License as published by
6 ** the Free Software Foundation; either version 2 of the License, or
7 ** (at your option) any later version.
8 **
9 ** This program is distributed in the hope that it will be useful,
10 ** but WITHOUT ANY WARRANTY; without even the implied warranty of
11 ** MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
12 ** GNU General Public License for more details.
13 **
14 ** You should have received a copy of the GNU General Public License
15 ** along with this program; if not, write to the Free Software
16 ** Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.
17 */
18
19 #include <stdio.h>
20 #include <stdlib.h>
21 #include <math.h>
22
23 #include "dft_cmp.h"
24 #include "utils.h"
25
26 #ifndef M_PI
27 #define M_PI 3.14159265358979323846264338
28 #endif
29
30 #define DFT_SPEC_LENGTH (DFT_DATA_LENGTH / 2)
31
32 static void dft_magnitude (const double *data, double *spectrum) ;
33 static double calc_max_spectral_difference (const double *spec1, const double *spec2) ;
34
35 /*--------------------------------------------------------------------------------
36 ** Public functions.
37 */
38
39 double
40 dft_cmp_float (int linenum, const float *in_data, const float *test_data, int len, double target_snr, int allow_exit)
41 { static double orig [DFT_DATA_LENGTH] ;
42 static double test [DFT_DATA_LENGTH] ;
43 unsigned k ;
44
45 if (len != DFT_DATA_LENGTH)
46 { printf ("Error (line %d) : dft_cmp_float : Bad input array length.\n", linenum) ;
47 return 1 ;
48 } ;
49
50 for (k = 0 ; k < ARRAY_LEN (orig) ; k++)
51 { test [k] = test_data [k] ;
52 orig [k] = in_data [k] ;
53 } ;
54
55 return dft_cmp_double (linenum, orig, test, len, target_snr, allow_exit) ;
56 } /* dft_cmp_float */
57
58 double
59 dft_cmp_double (int linenum, const double *orig, const double *test, int len, double target_snr, int allow_exit)
60 { static double orig_spec [DFT_SPEC_LENGTH] ;
61 static double test_spec [DFT_SPEC_LENGTH] ;
62 double snr ;
63
64 if (! orig || ! test)
65 { printf ("Error (line %d) : dft_cmp_double : Bad input arrays.\n", linenum) ;
66 return 1 ;
67 } ;
68
69 if (len != DFT_DATA_LENGTH)
70 { printf ("Error (line %d) : dft_cmp_double : Bad input array length.\n", linenum) ;
71 return 1 ;
72 } ;
73
74 dft_magnitude (orig, orig_spec) ;
75 dft_magnitude (test, test_spec) ;
76
77 snr = calc_max_spectral_difference (orig_spec, test_spec) ;
78
79 if (snr > target_snr)
80 { printf ("\n\nLine %d: Actual SNR (% 4.1f) > target SNR (% 4.1f).\n\n", linenum, snr, target_snr) ;
81 oct_save_double (orig, test, len) ;
82 if (allow_exit)
83 exit (1) ;
84 } ;
85
86 if (snr < -500.0)
87 snr = -500.0 ;
88
89 return snr ;
90 } /* dft_cmp_double */
91
92 /*--------------------------------------------------------------------------------
93 ** Quick dirty calculation of magnitude spectrum for real valued data using
94 ** Discrete Fourier Transform. Since the data is real, the DFT is only
95 ** calculated for positive frequencies.
96 */
97
98 static void
99 dft_magnitude (const double *data, double *spectrum)
100 { static double cos_angle [DFT_DATA_LENGTH] = { 0.0 } ;
101 static double sin_angle [DFT_DATA_LENGTH] ;
102
103 double real_part, imag_part ;
104 int k, n ;
105
106 /* If sine and cosine tables haven't been initialised, do so. */
107 if (cos_angle [0] == 0.0)
108 for (n = 0 ; n < DFT_DATA_LENGTH ; n++)
109 { cos_angle [n] = cos (2.0 * M_PI * n / DFT_DATA_LENGTH) ;
110 sin_angle [n] = -1.0 * sin (2.0 * M_PI * n / DFT_DATA_LENGTH) ;
111 } ;
112
113 /* DFT proper. Since the data is real, only generate a half spectrum. */
114 for (k = 1 ; k < DFT_SPEC_LENGTH ; k++)
115 { real_part = 0.0 ;
116 imag_part = 0.0 ;
117
118 for (n = 0 ; n < DFT_DATA_LENGTH ; n++)
119 { real_part += data [n] * cos_angle [(k * n) % DFT_DATA_LENGTH] ;
120 imag_part += data [n] * sin_angle [(k * n) % DFT_DATA_LENGTH] ;
121 } ;
122
123 spectrum [k] = sqrt (real_part * real_part + imag_part * imag_part) ;
124 } ;
125
126 spectrum [k] = 0.0 ;
127
128 spectrum [0] = spectrum [1] = spectrum [2] = spectrum [3] = spectrum [4] = 0.0 ;
129
130 return ;
131 } /* dft_magnitude */
132
133 static double
134 calc_max_spectral_difference (const double *orig, const double *test)
135 { double orig_max = 0.0, max_diff = 0.0 ;
136 int k ;
137
138 for (k = 0 ; k < DFT_SPEC_LENGTH ; k++)
139 { if (orig_max < orig [k])
140 orig_max = orig [k] ;
141 if (max_diff < fabs (orig [k] - test [k]))
142 max_diff = fabs (orig [k] - test [k]) ;
143 } ;
144
145 if (max_diff < 1e-25)
146 return -500.0 ;
147
148 return 20.0 * log10 (max_diff / orig_max) ;
149 } /* calc_max_spectral_difference */