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 };
11 #define NUM_COEFF (sizeof(coeffs)/sizeof(coeffs[0]))
12 #define FILTER_LAG (NUM_COEFF - 2)
17 if (fread(&x, 2, 1, in) != 1) {
23 short trunc_short(double x)
27 } else if (x <= -32768.0) {
34 void write_sample(double x)
36 short s = trunc_short(x * 0.5);
37 if (fwrite(&s, 2, 1, stdout) != 1) {
43 double filter_lo(double *s) {
45 for (unsigned i = 0; i < NUM_COEFF; ++i) {
46 acc += s[i] * coeffs[i];
51 double filter_hi(double *s) {
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];
60 void inv_filter(double *dst, double *src_lo, double *src_hi)
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];
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];
83 void filter_lo_series(double *dst, double *src, unsigned len)
85 for (unsigned i = 0; i < FILTER_LAG; ++i) {
86 dst[i] = dst[len + i];
88 for (unsigned i = 0; i < len; ++i) {
89 dst[i + FILTER_LAG] = filter_lo(src + i*2);
93 void filter_hi_series(double *dst, double *src, unsigned len)
95 for (unsigned i = 0; i < FILTER_LAG; ++i) {
96 dst[i] = dst[len + i];
98 for (unsigned i = 0; i < len; ++i) {
99 dst[i + FILTER_LAG] = filter_hi(src + i*2);
103 int main(int argc, char **argv)
105 in = fopen(argv[1], "rb");
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
119 for (unsigned i = 0; i < FILTER_LAG; ++i) {
120 s0[i] = s0[(1 << 4) + i];
122 for (unsigned i = 0; i < (1 << 4); ++i) {
123 s0[i + FILTER_LAG] = read_sample();
126 // First level (two bands)
127 filter_lo_series(s1[0], s0, (1 << 3));
128 filter_hi_series(s1[1], s0, (1 << 3));
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));
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));
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];
148 s4[i * 2][FILTER_LAG/2] = filter_lo(s3[i]);
149 s4[i * 2 + 1][FILTER_LAG/2] = filter_hi(s3[i]);
152 // COMPRESS HERE (based on s4[0..15][FILTER_LAG/2])
153 for (unsigned i = 0; i < 16; ++i) {
155 s4[i][FILTER_LAG/2] = 0.0f;
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];
163 inv_filter(u3[i] + FILTER_LAG/2, s4[i * 2], s4[i * 2 + 1]);
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];
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);
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];
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);
187 for (unsigned j = 0; j < 8; ++j) {
188 inv_filter(u0 + j*2, u1[0] + j, u1[1] + j);
191 for (unsigned j = 0; j < 16; ++j) {
198 for (unsigned j = 2; j < 18; ++j) {
199 printf("%f ", s0[j]);
202 for (unsigned j = 0; j < 16; ++j) {
203 printf("%f ", u0[j]);