c902769f817501a45b2bd9a90605d6b26ecf4958
2 ** Copyright (C) 2002-2011 Erik de Castro Lopo <erikd@mega-nerd.com>
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.
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.
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.
27 #define M_PI 3.14159265358979323846264338
30 #define DFT_SPEC_LENGTH (DFT_DATA_LENGTH / 2)
32 static void dft_magnitude (const double *data
, double *spectrum
) ;
33 static double calc_max_spectral_difference (const double *spec1
, const double *spec2
) ;
35 /*--------------------------------------------------------------------------------
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
] ;
45 if (len
!= DFT_DATA_LENGTH
)
46 { printf ("Error (line %d) : dft_cmp_float : Bad input array length.\n", linenum
) ;
50 for (k
= 0 ; k
< ARRAY_LEN (orig
) ; k
++)
51 { test
[k
] = test_data
[k
] ;
52 orig
[k
] = in_data
[k
] ;
55 return dft_cmp_double (linenum
, orig
, test
, len
, target_snr
, allow_exit
) ;
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
] ;
65 { printf ("Error (line %d) : dft_cmp_double : Bad input arrays.\n", linenum
) ;
69 if (len
!= DFT_DATA_LENGTH
)
70 { printf ("Error (line %d) : dft_cmp_double : Bad input array length.\n", linenum
) ;
74 dft_magnitude (orig
, orig_spec
) ;
75 dft_magnitude (test
, test_spec
) ;
77 snr
= calc_max_spectral_difference (orig_spec
, test_spec
) ;
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
) ;
90 } /* dft_cmp_double */
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.
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
] ;
103 double real_part
, imag_part
;
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
) ;
113 /* DFT proper. Since the data is real, only generate a half spectrum. */
114 for (k
= 1 ; k
< DFT_SPEC_LENGTH
; k
++)
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
] ;
123 spectrum
[k
] = sqrt (real_part
* real_part
+ imag_part
* imag_part
) ;
128 spectrum
[0] = spectrum
[1] = spectrum
[2] = spectrum
[3] = spectrum
[4] = 0.0 ;
131 } /* dft_magnitude */
134 calc_max_spectral_difference (const double *orig
, const double *test
)
135 { double orig_max
= 0.0, max_diff
= 0.0 ;
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
]) ;
145 if (max_diff
< 1e-25)
148 return 20.0 * log10 (max_diff
/ orig_max
) ;
149 } /* calc_max_spectral_difference */