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 */