--- /dev/null
+#include <stdio.h>
+#include <vector>
+
+FILE *in;
+
+//double coeffs[] = { 0.6830127, 1.1830127, 0.3169873, -0.1830127 };
+//double coeffs[] = { 0.47046721, 1.14111692, 0.650365, -0.19093442, -0.12083221, 0.0498175 };
+double coeffs[] = { 0.22641898, 0.85394354, 1.02432694, 0.19576696, -0.34265671, -0.04560113, 0.10970265, -0.00882680, -0.01779187, 4.71742793e-3 };
+//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 };
+
+#define NUM_COEFF (sizeof(coeffs)/sizeof(coeffs[0]))
+#define FILTER_LAG (NUM_COEFF - 2)
+
+double read_sample()
+{
+ short x;
+ if (fread(&x, 2, 1, in) != 1) {
+ exit(0);
+ }
+ return x;
+}
+
+short trunc_short(double x)
+{
+ if (x >= 32767.0) {
+ return 32767;
+ } else if (x <= -32768.0) {
+ return -32768;
+ } else {
+ return short(x);
+ }
+}
+
+void write_sample(double x)
+{
+ short s = trunc_short(x * 0.5);
+ if (fwrite(&s, 2, 1, stdout) != 1) {
+ perror("fwrite");
+ exit(1);
+ }
+}
+
+double filter_lo(double *s) {
+ double acc = 0.0;
+ for (unsigned i = 0; i < NUM_COEFF; ++i) {
+ acc += s[i] * coeffs[i];
+ }
+ return acc;
+}
+
+double filter_hi(double *s) {
+ double acc = 0.0;
+ for (unsigned i = 0; i < NUM_COEFF; i += 2) {
+ acc += s[i] * coeffs[NUM_COEFF - i - 1];
+ acc -= s[i + 1] * coeffs[NUM_COEFF - i - 2];
+ }
+ return acc;
+}
+
+void inv_filter(double *dst, double *src_lo, double *src_hi)
+{
+ // first sample
+ {
+ double c1 = 0.0, c2 = 0.0;
+ for (unsigned i = 0; i < NUM_COEFF/2; ++i) {
+ c1 += src_lo[i] * coeffs[NUM_COEFF - (i*2 + 1) - 1];
+ c2 += src_hi[i] * coeffs[i*2 + 1];
+ }
+ dst[0] = 0.5*(c1+c2);
+ }
+
+ // second sample
+ {
+ double c1 = 0.0, c2 = 0.0;
+ for (unsigned i = 0; i < NUM_COEFF/2; ++i) {
+ c1 += src_lo[i] * coeffs[NUM_COEFF - i*2 - 1];
+ c2 += src_hi[i] * coeffs[i*2];
+ }
+ dst[1] = 0.5*(c1-c2);
+ }
+}
+
+void filter_lo_series(double *dst, double *src, unsigned len)
+{
+ for (unsigned i = 0; i < FILTER_LAG; ++i) {
+ dst[i] = dst[len + i];
+ }
+ for (unsigned i = 0; i < len; ++i) {
+ dst[i + FILTER_LAG] = filter_lo(src + i*2);
+ }
+}
+
+void filter_hi_series(double *dst, double *src, unsigned len)
+{
+ for (unsigned i = 0; i < FILTER_LAG; ++i) {
+ dst[i] = dst[len + i];
+ }
+ for (unsigned i = 0; i < len; ++i) {
+ dst[i + FILTER_LAG] = filter_hi(src + i*2);
+ }
+}
+
+int main(int argc, char **argv)
+{
+ in = fopen(argv[1], "rb");
+
+ double s0[(1 << 4) + FILTER_LAG] = {0}; // 18
+ double s1[2][(1 << 3) + FILTER_LAG] = {{0}}; // 10
+ double s2[4][(1 << 2) + FILTER_LAG] = {{0}}; // 6
+ double s3[8][(1 << 1) + FILTER_LAG] = {{0}}; // 4
+ double s4[16][1 + FILTER_LAG/2] = {{0}}; // 2
+ double u3[8][(1 << 1) + FILTER_LAG/2] = {{0}}; // 3
+ double u2[4][(1 << 2) + FILTER_LAG/2] = {{0}}; // 5
+ double u1[2][(1 << 3) + FILTER_LAG/2] = {{0}}; // 9
+ double u0[1 << 4] = {0}; // 16
+
+ for ( ;; ) {
+ // Input (one band)
+ for (unsigned i = 0; i < FILTER_LAG; ++i) {
+ s0[i] = s0[(1 << 4) + i];
+ }
+ for (unsigned i = 0; i < (1 << 4); ++i) {
+ s0[i + FILTER_LAG] = read_sample();
+ }
+
+ // First level (two bands)
+ filter_lo_series(s1[0], s0, (1 << 3));
+ filter_hi_series(s1[1], s0, (1 << 3));
+
+ // Second level (four bands)
+ for (unsigned i = 0; i < 2; ++i) {
+ filter_lo_series(s2[i * 2], s1[i], (1 << 2));
+ filter_hi_series(s2[i * 2 + 1], s1[i], (1 << 2));
+ }
+
+ // Third level (eight bands)
+ for (unsigned i = 0; i < 4; ++i) {
+ filter_lo_series(s3[i * 2], s2[i], (1 << 1));
+ filter_hi_series(s3[i * 2 + 1], s2[i], (1 << 1));
+ }
+
+ // Fourth level (sixteen bands)
+ for (unsigned i = 0; i < 8; ++i) {
+ for (unsigned j = 0; j < FILTER_LAG/2; ++j) {
+ s4[i * 2][j] = s4[i * 2][j + 1];
+ s4[i * 2 + 1][j] = s4[i * 2 + 1][j + 1];
+ }
+ s4[i * 2][FILTER_LAG/2] = filter_lo(s3[i]);
+ s4[i * 2 + 1][FILTER_LAG/2] = filter_hi(s3[i]);
+ }
+
+ // COMPRESS HERE (based on s4[0..15][FILTER_LAG/2])
+ for (unsigned i = 0; i < 16; ++i) {
+ if (i != 0)
+ s4[i][FILTER_LAG/2] = 0.0f;
+ }
+
+ // Third level (eight bands)
+ for (unsigned i = 0; i < 8; ++i) {
+ for (unsigned j = 0; j < FILTER_LAG/2; ++j) {
+ u3[i][j] = u3[i][j + 2];
+ }
+ inv_filter(u3[i] + FILTER_LAG/2, s4[i * 2], s4[i * 2 + 1]);
+ }
+
+ // Second level (four bands)
+ for (unsigned i = 0; i < 4; ++i) {
+ for (unsigned j = 0; j < FILTER_LAG/2; ++j) {
+ u2[i][j] = u2[i][j + 4];
+ }
+ for (unsigned j = 0; j < 2; ++j) {
+ inv_filter(u2[i] + j*2 + FILTER_LAG/2, u3[i * 2] + j, u3[i * 2 + 1] + j);
+ }
+ }
+
+ // First level (two bands)
+ for (unsigned i = 0; i < 2; ++i) {
+ for (unsigned j = 0; j < FILTER_LAG/2; ++j) {
+ u1[i][j] = u1[i][j + 8];
+ }
+ for (unsigned j = 0; j < 4; ++j) {
+ inv_filter(u1[i] + j*2 + FILTER_LAG/2, u2[i * 2] + j, u2[i * 2 + 1] + j);
+ }
+ }
+
+ // Output data
+ for (unsigned j = 0; j < 8; ++j) {
+ inv_filter(u0 + j*2, u1[0] + j, u1[1] + j);
+ }
+
+ for (unsigned j = 0; j < 16; ++j) {
+ write_sample(u0[j]);
+ }
+
+#if 0
+ // DEBUG
+ printf("IN: ");
+ for (unsigned j = 2; j < 18; ++j) {
+ printf("%f ", s0[j]);
+ }
+ printf("\nOUT: ");
+ for (unsigned j = 0; j < 16; ++j) {
+ printf("%f ", u0[j]);
+ }
+ printf("\n\n");
+#endif
+ }
+}
projects / qmfsplit / commitdiff