freqdom.c (1299B)
1 #include <u.h> 2 #include <libc.h> 3 4 #include "freqdom.h" 5 6 #define dct_cos(X, U) dct_matrix[U * DCT_N + X] 7 8 s16int cos16_table[0x4000]; 9 10 double dct_alpha0, dct_alpha1, *dct_matrix; 11 12 void 13 init_cos16_table(void) 14 { 15 double i; 16 s16int *p = cos16_table; 17 for (i = 0; p < cos16_table + 0x4000; i += 1.0/(0x4000), p++) { 18 *p = (s16int)(cos(i * PIO2) * 0x7fff); 19 } 20 } 21 22 s16int 23 cos16(u16int t) 24 { 25 switch((t/0x4000)){ 26 default: 27 case 0: 28 return cos16_table[t&0x3fff]; 29 case 1: 30 return -cos16_table[0x3fff - t&0x3fff]; 31 case 2: 32 return -cos16_table[t&0x3fff]; 33 case 3: 34 return cos16_table[0x3fff - t&0x3fff]; 35 } 36 } 37 38 void 39 init_dct(void) 40 { 41 dct_matrix = malloc(sizeof(double) * DCT_N * DCT_N); 42 int x, u; 43 for (x = 0; x < DCT_N; x++) { 44 for (u = 0; u < DCT_N; u++) { 45 dct_cos(x, u) = cos(PI * ( 2 * x + 1) * u / (2 * DCT_N)); 46 } 47 } 48 dct_alpha0=sqrt(1.0/DCT_N); 49 dct_alpha1=sqrt(2.0/DCT_N); 50 } 51 52 s16int 53 dct(int u, s16int *buf) 54 { 55 int x; 56 double f, d = 0; 57 for (x = 0; x < DCT_N; x++) { 58 f = buf[x]; 59 d += f * dct_cos(x, u); 60 } 61 if (u == 0) return dct_alpha0 * d; 62 return dct_alpha1 * d; 63 } 64 65 s16int 66 invdct(int x, s16int *buf) 67 { 68 int u; 69 double f, p, d = 0; 70 double alpha = dct_alpha0; 71 for (u = 0; u < DCT_N; u++) { 72 f = buf[u]; 73 p = alpha * f * dct_cos(x, u); 74 d += p; 75 alpha = dct_alpha1; 76 } 77 return d; 78 }