git.sesse.net Git - qmfsplit/blob - qmfsplit.cpp

   1 #include <stdio.h>
   2 #include <vector>
   3
   4 FILE *in;
   5
   6 //double coeffs[] = { 0.6830127, 1.1830127, 0.3169873, -0.1830127 };
   7 //double coeffs[] = { 0.47046721, 1.14111692, 0.650365, -0.19093442, -0.12083221, 0.0498175 };
   8 double coeffs[] = { 0.22641898, 0.85394354, 1.02432694, 0.19576696, -0.34265671, -0.04560113, 0.10970265, -0.00882680, -0.01779187, 4.71742793e-3 };
   9 //double coeffs[] = { 0.03771716, 0.26612218, 0.74557507, 0.97362811, 0.39763774, -0.35333620, -0.27710988, 0.18012745, 0.13160299, -0.10096657, -0.04165925, 0.04696981, 5.10043697e-3, -0.01517900, 1.97332536e-3, 2.81768659e-3, -9.69947840e-4, -1.64709006e-4, 1.32354367e-4, -1.875841e-5 };
  10
  11 #define NUM_COEFF (sizeof(coeffs)/sizeof(coeffs[0]))
  12 #define FILTER_LAG (NUM_COEFF - 2)
  13
  14 double read_sample()
  15 {
  16         short x;
  17         if (fread(&x, 2, 1, in) != 1) {
  18                 exit(0);
  19         }
  20         return x;
  21 }
  22
  23 short trunc_short(double x)
  24 {
  25         if (x >= 32767.0) {
  26                 return 32767;
  27         } else if (x <= -32768.0) {
  28                 return -32768;
  29         } else {
  30                 return short(x);
  31         }
  32 }
  33
  34 void write_sample(double x)
  35 {
  36         short s = trunc_short(x * 0.5);
  37         if (fwrite(&s, 2, 1, stdout) != 1) {
  38                 perror("fwrite");
  39                 exit(1);
  40         }
  41 }
  42
  43 double filter_lo(double *s) {
  44         double acc = 0.0;
  45         for (unsigned i = 0; i < NUM_COEFF; ++i) {
  46                 acc += s[i] * coeffs[i];
  47         }
  48         return acc;
  49 }
  50
  51 double filter_hi(double *s) {
  52         double acc = 0.0;
  53         for (unsigned i = 0; i < NUM_COEFF; i += 2) {
  54                 acc += s[i] * coeffs[NUM_COEFF - i - 1];
  55                 acc -= s[i + 1] * coeffs[NUM_COEFF - i - 2];
  56         }
  57         return acc;
  58 }
  59
  60 void inv_filter(double *dst, double *src_lo, double *src_hi)
  61 {
  62         // first sample
  63         {
  64                 double c1 = 0.0, c2 = 0.0;
  65                 for (unsigned i = 0; i < NUM_COEFF/2; ++i) {
  66                         c1 += src_lo[i] * coeffs[NUM_COEFF - (i*2 + 1) - 1];
  67                         c2 += src_hi[i] * coeffs[i*2 + 1];
  68                 }
  69                 dst[0] = 0.5*(c1+c2);
  70         }
  71
  72         // second sample
  73         {
  74                 double c1 = 0.0, c2 = 0.0;
  75                 for (unsigned i = 0; i < NUM_COEFF/2; ++i) {
  76                         c1 += src_lo[i] * coeffs[NUM_COEFF - i*2 - 1];
  77                         c2 += src_hi[i] * coeffs[i*2];
  78                 }
  79                 dst[1] = 0.5*(c1-c2);
  80         }
  81 }
  82
  83 void filter_lo_series(double *dst, double *src, unsigned len)
  84 {
  85         for (unsigned i = 0; i < FILTER_LAG; ++i) {
  86                 dst[i] = dst[len + i];
  87         }
  88         for (unsigned i = 0; i < len; ++i) {
  89                 dst[i + FILTER_LAG] = filter_lo(src + i*2);
  90         }
  91 }
  92
  93 void filter_hi_series(double *dst, double *src, unsigned len)
  94 {
  95         for (unsigned i = 0; i < FILTER_LAG; ++i) {
  96                 dst[i] = dst[len + i];
  97         }
  98         for (unsigned i = 0; i < len; ++i) {
  99                 dst[i + FILTER_LAG] = filter_hi(src + i*2);
 100         }
 101 }
 102
 103 int main(int argc, char **argv)
 104 {
 105         in = fopen(argv[1], "rb");
 106
 107         double s0[(1 << 4) + FILTER_LAG] = {0};           // 18
 108         double s1[2][(1 << 3) + FILTER_LAG] = {{0}};      // 10
 109         double s2[4][(1 << 2) + FILTER_LAG] = {{0}};      // 6
 110         double s3[8][(1 << 1) + FILTER_LAG] = {{0}};      // 4
 111         double s4[16][1 + FILTER_LAG/2] = {{0}};          // 2
 112         double u3[8][(1 << 1) + FILTER_LAG/2] = {{0}};    // 3
 113         double u2[4][(1 << 2) + FILTER_LAG/2] = {{0}};    // 5
 114         double u1[2][(1 << 3) + FILTER_LAG/2] = {{0}};    // 9
 115         double u0[1 << 4] = {0};                          // 16
 116
 117         for ( ;; ) {
 118                 // Input (one band)
 119                 for (unsigned i = 0; i < FILTER_LAG; ++i) {
 120                         s0[i] = s0[(1 << 4) + i];
 121                 }
 122                 for (unsigned i = 0; i < (1 << 4); ++i) {
 123                         s0[i + FILTER_LAG] = read_sample();
 124                 }
 125
 126                 // First level (two bands)
 127                 filter_lo_series(s1[0], s0, (1 << 3));
 128                 filter_hi_series(s1[1], s0, (1 << 3));
 129
 130                 // Second level (four bands)
 131                 for (unsigned i = 0; i < 2; ++i) {
 132                         filter_lo_series(s2[i * 2],     s1[i], (1 << 2));
 133                         filter_hi_series(s2[i * 2 + 1], s1[i], (1 << 2));
 134                 }
 135
 136                 // Third level (eight bands)
 137                 for (unsigned i = 0; i < 4; ++i) {
 138                         filter_lo_series(s3[i * 2],     s2[i], (1 << 1));
 139                         filter_hi_series(s3[i * 2 + 1], s2[i], (1 << 1));
 140                 }
 141
 142                 // Fourth level (sixteen bands)
 143                 for (unsigned i = 0; i < 8; ++i) {
 144                         for (unsigned j = 0; j < FILTER_LAG/2; ++j) {
 145                                 s4[i * 2][j]     = s4[i * 2][j + 1];
 146                                 s4[i * 2 + 1][j] = s4[i * 2 + 1][j + 1];
 147                         }
 148                         s4[i * 2][FILTER_LAG/2]     = filter_lo(s3[i]);
 149                         s4[i * 2 + 1][FILTER_LAG/2] = filter_hi(s3[i]);
 150                 }
 151
 152                 // COMPRESS HERE (based on s4[0..15][FILTER_LAG/2])
 153                 for (unsigned i = 0; i < 16; ++i) {
 154                         if (i != 0)
 155                                 s4[i][FILTER_LAG/2] = 0.0f;
 156                 }
 157
 158                 // Third level (eight bands)
 159                 for (unsigned i = 0; i < 8; ++i) {
 160                         for (unsigned j = 0; j < FILTER_LAG/2; ++j) {
 161                                 u3[i][j] = u3[i][j + 2];
 162                         }
 163                         inv_filter(u3[i] + FILTER_LAG/2, s4[i * 2], s4[i * 2 + 1]);
 164                 }
 165
 166                 // Second level (four bands)
 167                 for (unsigned i = 0; i < 4; ++i) {
 168                         for (unsigned j = 0; j < FILTER_LAG/2; ++j) {
 169                                 u2[i][j] = u2[i][j + 4];
 170                         }
 171                         for (unsigned j = 0; j < 2; ++j) {
 172                                 inv_filter(u2[i] + j*2 + FILTER_LAG/2, u3[i * 2] + j, u3[i * 2 + 1] + j);
 173                         }
 174                 }
 175
 176                 // First level (two bands)
 177                 for (unsigned i = 0; i < 2; ++i) {
 178                         for (unsigned j = 0; j < FILTER_LAG/2; ++j) {
 179                                 u1[i][j] = u1[i][j + 8];
 180                         }
 181                         for (unsigned j = 0; j < 4; ++j) {
 182                                 inv_filter(u1[i] + j*2 + FILTER_LAG/2, u2[i * 2] + j, u2[i * 2 + 1] + j);
 183                         }
 184                 }
 185
 186                 // Output data
 187                 for (unsigned j = 0; j < 8; ++j) {
 188                         inv_filter(u0 + j*2, u1[0] + j, u1[1] + j);
 189                 }
 190
 191                 for (unsigned j = 0; j < 16; ++j) {
 192                         write_sample(u0[j]);
 193                 }
 194
 195 #if 0
 196                 // DEBUG
 197                 printf("IN: ");
 198                 for (unsigned j = 2; j < 18; ++j) {
 199                         printf("%f ", s0[j]);
 200                 }
 201                 printf("\nOUT: ");
 202                 for (unsigned j = 0; j < 16; ++j) {
 203                         printf("%f ", u0[j]);
 204                 }
 205                 printf("\n\n");
 206 #endif
 207         }
 208 }