testsuite/wavelet_noise_gain.c

   1
   2 #ifdef HAVE_CONFIG_H
   3 #include "config.h"
   4 #endif
   5
   6 #include <schroedinger/schro.h>
   7 #include <schroedinger/schrowavelet.h>
   8 #include <schroedinger/schrooil.h>
   9 #include <liboil/liboil.h>
  10
  11 #include <stdio.h>
  12 #include <string.h>
  13 #include <stdlib.h>
  14 #include <math.h>
  15
  16 #include "common.h"
  17
  18 int16_t tmp[2000];
  19 int16_t tmp2[2000];
  20
  21 int filtershift[] = { 1, 1, 1, 0, 1, 0, 1 };
  22
  23 void synth(int16_t *a, int filter, int n);
  24 void split (int16_t *a, int filter, int n);
  25 void synth_schro_ext (int16_t *a, int filter, int n);
  26 void split_schro_ext (int16_t *a, int n, int filter);
  27 void deinterleave (int16_t *a, int n);
  28 void interleave (int16_t *a, int n);
  29 void dump (int16_t *a, int n);
  30 void dump_cmp (int16_t *a, int16_t *b, int n);
  31
  32 #define SHIFT 8
  33 #define N (1<<SHIFT)
  34
  35 double
  36 gain_test(int filter, int level, int n_levels)
  37 {
  38   int16_t *a = tmp+10;
  39   double gain;
  40   int i;
  41   int j;
  42   double p1,p2;
  43
  44   gain = 0;
  45
  46   for(j=0;j<10000;j++){
  47     p1 = 0;
  48     p2 = 0;
  49     for(i=0;i<N;i++){
  50       a[i] = 0;
  51     }
  52     if (level == 0) {
  53       for(i=0;i<N>>n_levels;i++){
  54         a[i]=rint(100 * random_std());
  55       }
  56     } else {
  57       for(i=N>>(n_levels-level+1);i<N>>(n_levels-level);i++){
  58         a[i]=rint(100 * random_std());
  59       }
  60     }
  61     for(i=0;i<N;i++){
  62       p1 += a[i]*a[i];
  63     }
  64
  65     for(i=0;i<n_levels;i++){
  66       interleave(a,N>>(n_levels-1-i));
  67       synth_schro_ext(a,N>>(n_levels-1-i),filter);
  68     }
  69
  70     for(i=0;i<256;i++){
  71       p2 += a[i]*a[i];
  72     }
  73
  74     gain += p2/p1;
  75   }
  76
  77   //printf("%d %g %g\n", level, gain/10000, sqrt(gain/10000));
  78
  79   return sqrt(gain/10000);
  80 }
  81
  82 int
  83 main (int argc, char *argv[])
  84 {
  85   int filter;
  86   int n_levels;
  87
  88   schro_init();
  89
  90   filter = 0;
  91   if (argc > 1) {
  92     filter = strtol(argv[1], NULL, 0);
  93   }
  94
  95   {
  96     double a[5];
  97     double b[9];
  98     int i;
  99     double min;
 100
 101     for(n_levels=1;n_levels<=4;n_levels++){
 102       printf("n_levels=%d:\n", n_levels);
 103       for(i=0;i<=n_levels;i++){
 104         a[i] = gain_test (filter, i, n_levels);
 105       }
 106       for(i=0;i<=n_levels;i++){
 107         printf("%d %5.3f\n", i, a[i]);
 108       }
 109
 110       b[0] = a[0] * a[0] / (1<<(n_levels*filtershift[filter]));
 111       for(i=0;i<n_levels;i++){
 112         b[i*2+1] = a[i+0] * a[i+1] / (1<<((n_levels-i)*filtershift[filter]));
 113         b[i*2+2] = a[i+1] * a[i+1] / (1<<((n_levels-i)*filtershift[filter]));
 114       }
 115
 116       min = b[0];
 117       for(i=0;i<n_levels*2+1;i++){
 118         if (b[i] < min) min = b[i];
 119       }
 120
 121       for(i=0;i<n_levels*2+1;i++){
 122         printf("%d %5.3f %5.3f %d\n", i, b[i], b[i]/min, gain_to_quant_index(b[i]/min));
 123       }
 124     }
 125   }
 126
 127   return 0;
 128 }
 129
 130 void
 131 interleave (int16_t *a, int n)
 132 {
 133   int i;
 134   for(i=0;i<n/2;i++){
 135     tmp2[i*2] = a[i];
 136     tmp2[i*2 + 1] = a[n/2 + i];
 137   }
 138   for(i=0;i<n;i++){
 139     a[i] = tmp2[i];
 140   }
 141 }
 142
 143 void
 144 deinterleave (int16_t *a, int n)
 145 {
 146   int i;
 147   for(i=0;i<n/2;i++){
 148     tmp2[i] = a[i*2];
 149     tmp2[n/2 + i] = a[i*2+1];
 150   }
 151   for(i=0;i<n;i++){
 152     a[i] = tmp2[i];
 153   }
 154 }
 155
 156 void
 157 split_schro_ext (int16_t *a, int n, int filter)
 158 {
 159   int16_t tmp1[2000], *hi;
 160   int16_t tmp2[2000], *lo;
 161
 162   hi = tmp1 + 4;
 163   lo = tmp2 + 4;
 164
 165   oil_deinterleave2_s16 (hi, lo, a, n/2);
 166
 167   switch (filter) {
 168     case SCHRO_WAVELET_DESLAURIERS_DUBUC_9_7:
 169       schro_split_ext_desl93 (hi, lo, n/2);
 170       break;
 171     case SCHRO_WAVELET_LE_GALL_5_3:
 172       schro_split_ext_53 (hi, lo, n/2);
 173       break;
 174     case SCHRO_WAVELET_DESLAURIERS_DUBUC_13_7:
 175       schro_split_ext_135 (hi, lo, n/2);
 176       break;
 177     case SCHRO_WAVELET_HAAR_0:
 178     case SCHRO_WAVELET_HAAR_1:
 179       schro_split_ext_haar (hi, lo, n/2);
 180       break;
 181     case SCHRO_WAVELET_FIDELITY:
 182       schro_split_ext_fidelity (hi, lo, n/2);
 183       break;
 184     case SCHRO_WAVELET_DAUBECHIES_9_7:
 185       schro_split_ext_daub97(hi, lo, n/2);
 186       break;
 187   }
 188   oil_interleave2_s16 (a, hi, lo, n/2);
 189
 190 }
 191
 192 void
 193 synth_schro_ext (int16_t *a, int n, int filter)
 194 {
 195   int16_t tmp1[2000], *hi;
 196   int16_t tmp2[2000], *lo;
 197
 198   hi = tmp1 + 4;
 199   lo = tmp2 + 4;
 200
 201   oil_deinterleave2_s16 (hi, lo, a, n/2);
 202
 203   switch (filter) {
 204     case SCHRO_WAVELET_DESLAURIERS_DUBUC_9_7:
 205       schro_synth_ext_desl93 (hi, lo, n/2);
 206       break;
 207     case SCHRO_WAVELET_LE_GALL_5_3:
 208       schro_synth_ext_53 (hi, lo, n/2);
 209       break;
 210     case SCHRO_WAVELET_DESLAURIERS_DUBUC_13_7:
 211       schro_synth_ext_135 (hi, lo, n/2);
 212       break;
 213     case SCHRO_WAVELET_HAAR_0:
 214     case SCHRO_WAVELET_HAAR_1:
 215       schro_synth_ext_haar (hi, lo, n/2);
 216       break;
 217     case SCHRO_WAVELET_FIDELITY:
 218       schro_synth_ext_fidelity (hi, lo, n/2);
 219       break;
 220     case SCHRO_WAVELET_DAUBECHIES_9_7:
 221       schro_synth_ext_daub97(hi, lo, n/2);
 222       break;
 223   }
 224
 225   oil_interleave2_s16 (a, hi, lo, n/2);
 226 }
 227
 228
 229
 230