1 /*
2 ** Copyright (C) 2002-2015 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 "sfconfig.h"
20
21 #include <stdio.h>
22 #include <stdlib.h>
23 #include <math.h>
24
25 #include "dft_cmp.h"
26 #include "utils.h"
27
28 #ifndef M_PI
29 #define M_PI 3.14159265358979323846264338
30 #endif
31
32 #define DFT_SPEC_LENGTH (DFT_DATA_LENGTH / 2)
33
34 static void dft_magnitude (const double *data, double *spectrum) ;
35 static double calc_max_spectral_difference (const double *spec1, const double *spec2) ;
36
37 /*--------------------------------------------------------------------------------
38 ** Public functions.
39 */
40
41 double
dft_cmp_float(int linenum,const float * in_data,const float * test_data,int len,double target_snr,int allow_exit)42 dft_cmp_float (int linenum, const float *in_data, const float *test_data, int len, double target_snr, int allow_exit)
43 { static double orig [DFT_DATA_LENGTH] ;
44 static double test [DFT_DATA_LENGTH] ;
45 unsigned k ;
46
47 if (len != DFT_DATA_LENGTH)
48 { printf ("Error (line %d) : dft_cmp_float : Bad input array length.\n", linenum) ;
49 return 1 ;
50 } ;
51
52 for (k = 0 ; k < ARRAY_LEN (orig) ; k++)
53 { test [k] = test_data [k] ;
54 orig [k] = in_data [k] ;
55 } ;
56
57 return dft_cmp_double (linenum, orig, test, len, target_snr, allow_exit) ;
58 } /* dft_cmp_float */
59
60 double
dft_cmp_double(int linenum,const double * orig,const double * test,int len,double target_snr,int allow_exit)61 dft_cmp_double (int linenum, const double *orig, const double *test, int len, double target_snr, int allow_exit)
62 { static double orig_spec [DFT_SPEC_LENGTH] ;
63 static double test_spec [DFT_SPEC_LENGTH] ;
64 double snr ;
65
66 if (! orig || ! test)
67 { printf ("Error (line %d) : dft_cmp_double : Bad input arrays.\n", linenum) ;
68 return 1 ;
69 } ;
70
71 if (len != DFT_DATA_LENGTH)
72 { printf ("Error (line %d) : dft_cmp_double : Bad input array length.\n", linenum) ;
73 return 1 ;
74 } ;
75
76 dft_magnitude (orig, orig_spec) ;
77 dft_magnitude (test, test_spec) ;
78
79 snr = calc_max_spectral_difference (orig_spec, test_spec) ;
80
81 if (snr > target_snr)
82 { printf ("\n\nLine %d: Actual SNR (% 4.1f) > target SNR (% 4.1f).\n\n", linenum, snr, target_snr) ;
83 oct_save_double (orig, test, len) ;
84 if (allow_exit)
85 exit (1) ;
86 } ;
87
88 if (snr < -500.0)
89 snr = -500.0 ;
90
91 return snr ;
92 } /* dft_cmp_double */
93
94 /*--------------------------------------------------------------------------------
95 ** Quick dirty calculation of magnitude spectrum for real valued data using
96 ** Discrete Fourier Transform. Since the data is real, the DFT is only
97 ** calculated for positive frequencies.
98 */
99
100 static void
dft_magnitude(const double * data,double * spectrum)101 dft_magnitude (const double *data, double *spectrum)
102 { static double cos_angle [DFT_DATA_LENGTH] = { 0.0 } ;
103 static double sin_angle [DFT_DATA_LENGTH] ;
104
105 double real_part, imag_part ;
106 int k, n ;
107
108 /* If sine and cosine tables haven't been initialised, do so. */
109 if (cos_angle [0] == 0.0)
110 for (n = 0 ; n < DFT_DATA_LENGTH ; n++)
111 { cos_angle [n] = cos (2.0 * M_PI * n / DFT_DATA_LENGTH) ;
112 sin_angle [n] = -1.0 * sin (2.0 * M_PI * n / DFT_DATA_LENGTH) ;
113 } ;
114
115 /* DFT proper. Since the data is real, only generate a half spectrum. */
116 for (k = 1 ; k < DFT_SPEC_LENGTH ; k++)
117 { real_part = 0.0 ;
118 imag_part = 0.0 ;
119
120 for (n = 0 ; n < DFT_DATA_LENGTH ; n++)
121 { real_part += data [n] * cos_angle [(k * n) % DFT_DATA_LENGTH] ;
122 imag_part += data [n] * sin_angle [(k * n) % DFT_DATA_LENGTH] ;
123 } ;
124
125 spectrum [k] = sqrt (real_part * real_part + imag_part * imag_part) ;
126 } ;
127
128 spectrum [DFT_SPEC_LENGTH - 1] = 0.0 ;
129
130 spectrum [0] = spectrum [1] = spectrum [2] = 0.0 ;
131
132 return ;
133 } /* dft_magnitude */
134
135 static double
calc_max_spectral_difference(const double * orig,const double * test)136 calc_max_spectral_difference (const double *orig, const double *test)
137 { double orig_max = 0.0, max_diff = 0.0 ;
138 int k ;
139
140 for (k = 0 ; k < DFT_SPEC_LENGTH ; k++)
141 { if (orig_max < orig [k])
142 orig_max = orig [k] ;
143 if (max_diff < fabs (orig [k] - test [k]))
144 max_diff = fabs (orig [k] - test [k]) ;
145 } ;
146
147 if (max_diff < 1e-25)
148 return -500.0 ;
149
150 return 20.0 * log10 (max_diff / orig_max) ;
151 } /* calc_max_spectral_difference */
152