数据压缩之Polar编码实现详解

    Polar编码与范式Huffman编码类似，也是根据一个静态频率统计表来为各个符号分配前缀码来实现压缩。但Polar编码的构造算法不需要用Huffman的树结构，实现起来非常简单，而且大多数情况下压缩效果不会比Huffman差多少。
    Polar编码和范式Huffman编码的差别只在计算每个符号的编码长度上，这里举一个例子说明Polar编码的算法：

    假设拿到这样一个频率表：

    符号    频率
    A       190
    B       38
    C       185
    D       70
    E       253
    Total   736

    首先要按每个符号的频率降序排序：

    符号    频率
    E       253
    A       190
    C       185
    D       70
    B       38
    Total   736

    第二步，每个符号频率值按2的整数次幂向下取整，Total按2的整数次幂向上取整:

    符号    频率
    E       128
    A       128
    C       128
    D       64
    B       32
    Total   480
    Total‘  1024

    第三步也就是最关键的一步，进行前两步操作后，实际的频率总和Total必然不大于取整后的频率和Total'，现在我们从上到下看，对每个符号的频率，如果它倍增后实际总和不超过Total'，就对它进行倍增操作：

    符号    频率    倍增后
    E       128     256
    A       128     256
    C       128     256
    D       64      128
    B       32      64
    Total   480     960
    Total‘  1024

    第一轮操作，每个符号都满足了倍增的条件，而最终的实际频率和仍然不大于Total'，我们再进行一轮相同的操作：

    符号    频率    倍增后  第二轮倍增操作后
    E       128     256     256
    A       128     256     256
    C       128     256     256
    D       64      128     128
    B       32      64      128
    Total   480     960     1024
    Total‘  1024

    第二轮操作中只有最后一个符号满足倍增的条件，并且最后的实际频率和正好等于Total'，这时候操作结束，最后可以简单地用以下公式得出每个符号的编码长度：

  CodeLen(x) = Log2(Total / Freq(x))

用公式算出每个符号的编码长度：

    符号    频率    倍增后  第二轮倍增操作后    所求编码长度
    E       128     256     256                 Log2(1024 / 256) = 2
    A       128     256     256                 Log2(1024 / 256) = 2
    C       128     256     256                 Log2(1024 / 256) = 2
    D       64      128     128                 Log2(1024 / 128) = 3
    B       32      64      128                 Log2(1024 / 128) = 3
    Total   480     960     1024
    Total‘  1024

    至此Polar编码算法结束。

　有了每个符号的编码长度之后，我们就可以用范式Huffman编码类似的方法为每个符号分配前缀编码，进一步实现对数据的压缩了。

    最后提供一个完整的运用Polar编码进行无损数据压缩的简单实现，为了提高压缩率，我们对输入的数据先进行了一次MTF变换处理，使得最终的压缩率接近gzip。

  1 /*******************************************************************************
  2  * RichSelian's nooblike compressor: MTF + Polar coding
  3  ******************************************************************************/
  4 #include <stdio.h>
  5 #include <string.h>
  6 
  7 /*******************************************************************************
  8  * POLAR Coder
  9  ******************************************************************************/
 10 #define POLAR_SYMBOLS   256 /* should be even */
 11 #define POLAR_MAXLEN    15  /* should be less than 16, so we can pack two length values into a byte */
 12 
 13 #define M_round_down(x)     while((x)&(-(x)^(x))) { (x) &= (-(x)^(x)); }
 14 #define M_round_up(x)       while((x)&(-(x)^(x))) { (x) &= (-(x)^(x)); } (x) <<= 1;
 15 #define M_int_swap(x, y)    {int (_)=(x); (x)=(y); (y)=(_);}
 16 
 17 int polar_make_leng_table(const int* freq_table, int* leng_table) {
 18     int symbols[POLAR_SYMBOLS];
 19     int i;
 20     int s;
 21     int total;
 22     int shift = 0;
 23 
 24     memcpy(leng_table, freq_table, POLAR_SYMBOLS * sizeof(int));
 25 
 26 MakeTablePass:
 27     /* sort symbols */
 28     for(i = 0; i < POLAR_SYMBOLS; i++) {
 29         symbols[i] = i;
 30     }
 31     for(i = 0; i < POLAR_SYMBOLS; i++) {
 32         if(i > 0 && leng_table[symbols[i - 1]] < leng_table[symbols[i]]) {
 33             M_int_swap(symbols[i - 1], symbols[i]);
 34             i -= 2;
 35         }
 36     }
 37 
 38     /* calculate total frequency */
 39     total = 0;
 40     for(i = 0; i < POLAR_SYMBOLS; i++) {
 41         total += leng_table[i];
 42     }
 43 
 44     /* run */
 45     M_round_up(total);
 46     s = 0;
 47     for(i = 0; i < POLAR_SYMBOLS; i++) {
 48         M_round_down(leng_table[i]);
 49         s += leng_table[i];
 50     }
 51     while(s < total) {
 52         for(i = 0; i < POLAR_SYMBOLS; i++) {
 53             if(s + leng_table[symbols[i]] <= total) {
 54                 s += leng_table[symbols[i]];
 55                 leng_table[symbols[i]] *= 2;
 56             }
 57         }
 58     }
 59 
 60     /* get code length */
 61     for(i = 0; i < POLAR_SYMBOLS; i++) {
 62         s = 2;
 63         if(leng_table[i] > 0) {
 64             while((total / leng_table[i]) >> s != 0) {
 65                 s += 1;
 66             }
 67             leng_table[i] = s - 1;
 68         } else {
 69             leng_table[i] = 0;
 70         }
 71 
 72         /* code length too long -- scale and rebuild table */
 73         if(leng_table[i] > POLAR_MAXLEN) {
 74             shift += 1;
 75             for(i = 0; i < POLAR_SYMBOLS; i++) {
 76                 if((leng_table[i] = freq_table[i] >> shift) == 0 && freq_table[i] > 0) {
 77                     leng_table[i] = 1;
 78                 }
 79             }
 80             goto MakeTablePass;
 81         }
 82     }
 83     return 0;
 84 }
 85 
 86 int polar_make_code_table(const int* leng_table, int* code_table) {
 87     int i;
 88     int s;
 89     int t1;
 90     int t2;
 91     int code = 0;
 92 
 93     memset(code_table, 0, POLAR_SYMBOLS * sizeof(int));
 94 
 95     /* make code for each symbol */
 96     for(s = 1; s <= POLAR_MAXLEN; s++) {
 97         for(i = 0; i < POLAR_SYMBOLS; i++) {
 98             if(leng_table[i] == s) {
 99                 code_table[i] = code;
100                 code += 1;
101             }
102         }
103         code *= 2;
104     }
105 
106     /* reverse each code */
107     for(i = 0; i < POLAR_SYMBOLS; i++) {
108         t1 = 0;
109         t2 = leng_table[i] - 1;
110         while(t1 < t2) {
111             code_table[i] ^= (1 & (code_table[i] >> t1)) << t2;
112             code_table[i] ^= (1 & (code_table[i] >> t2)) << t1;
113             code_table[i] ^= (1 & (code_table[i] >> t1)) << t2;
114             t1++;
115             t2--;
116         }
117     }
118     return 0;
119 }
120 
121 int polar_make_decode_table(const int* leng_table, const int* code_table, int* decode_table) {
122     int i;
123     int c;
124 
125     for(c = 0; c < POLAR_SYMBOLS; c++) {
126         if(leng_table[c] > 0) {
127             for(i = 0; i + code_table[c] < 65536; i += (1 << leng_table[c])) {
128                 decode_table[i + code_table[c]] = c;
129             }
130         }
131     }
132     return 0;
133 }
134 
135 /*******************************************************************************
136  * MTF Transformer
137  ******************************************************************************/
138 unsigned char mtf_queue[65536][256];
139 unsigned char mtf_ch1 = 0;
140 unsigned char mtf_ch2 = 0;
141 unsigned char mtf_ch3 = 0;
142 
143 #define M_mtf_hashctx ((unsigned)(mtf_ch1*131313131 + mtf_ch2*13131 + mtf_ch3) % 65536)
144 
145 int mtf_transformer_init() {
146     int i;
147 
148     for(i = 0; i < sizeof(mtf_queue); i++) {
149         mtf_queue[i / 256][i % 256] = i % 256;
150     }
151     return 0;
152 }
153 
154 int mtf_encode(int c) {
155     int i;
156 
157     for(i = 0; i < 256; i++) {
158         if(mtf_queue[M_mtf_hashctx][i] == c) {
159             memmove(mtf_queue[M_mtf_hashctx] + 1, mtf_queue[M_mtf_hashctx], i);
160             mtf_queue[M_mtf_hashctx][0] = c;
161             break;
162         }
163     }
164     mtf_ch3 = mtf_ch2;
165     mtf_ch2 = mtf_ch1;
166     mtf_ch1 = c;
167     return i;
168 }
169 
170 int mtf_decode(int i) {
171     int c;
172 
173     c = mtf_queue[M_mtf_hashctx][i];
174     memmove(mtf_queue[M_mtf_hashctx] + 1, mtf_queue[M_mtf_hashctx], i);
175     mtf_queue[M_mtf_hashctx][0] = c;
176 
177     mtf_ch3 = mtf_ch2;
178     mtf_ch2 = mtf_ch1;
179     mtf_ch1 = c;
180     return c;
181 }
182 
183 /*******************************************************************************
184  * MAIN
185  ******************************************************************************/
186 int main(int argc, char** argv) {
187     static unsigned char idata[100000];
188     static unsigned char odata[120000]; /* never overflow */
189     int inpos;
190     int inlen;
191     int outpos;
192     int outlen;
193     int i;
194     int freq_table[256];
195     int code_table[256];
196     int leng_table[256];
197     int code_len = 0;
198     int code_buf = 0;
199     int decode_table[65536];
200 
201     mtf_transformer_init();
202 
203     /* compress */
204     if(argc == 2 && strcmp(argv[1], "e") == 0) {
205         /* init frequency table */
206         while((inlen = fread(idata, 1, sizeof(idata), stdin)) > 0) {
207             outlen = 0;
208             memset(freq_table, 0, sizeof(freq_table));
209 
210             for(i = 0; i < inlen; i++) {
211                 idata[i] = mtf_encode(idata[i]);
212                 freq_table[idata[i]] += 1;
213             }
214 
215             /* write length table */
216             polar_make_leng_table(freq_table, leng_table);
217             for(i = 0; i < POLAR_SYMBOLS; i += 2) {
218                 odata[outlen++] = leng_table[i] * 16 + leng_table[i + 1];
219             }
220 
221             /* encode */
222             polar_make_code_table(leng_table, code_table);
223             for(i = 0; i < inlen; i++) {
224                 code_buf += code_table[idata[i]] << code_len;
225                 code_len += leng_table[idata[i]];
226                 while(code_len > 8) {
227                     odata[outlen++] = code_buf % 256;
228                     code_buf /= 256;
229                     code_len -= 8;
230                 }
231             }
232             if(code_len > 0) {
233                 odata[outlen++] = code_buf;
234                 code_buf = 0;
235                 code_len = 0;
236             }
237             fwrite(&outlen, sizeof(outlen), 1, stdout);
238             fwrite(&inlen,  sizeof(inlen),  1, stdout);
239             fwrite(odata, 1, outlen, stdout);
240         }
241         return 0;
242     }
243 
244     /* decompress */
245     if(argc == 2 && strcmp(argv[1], "d") == 0) {
246         while(1) {
247             code_buf = 0;
248             code_len = 0;
249 
250             /* read inlen/outpos */
251             if(fread(&outlen, sizeof(outlen), 1, stdin) != 1 || fread(&inlen, sizeof(inlen), 1, stdin) != 1) {
252                 break;
253             }
254             fread(odata, 1, outlen, stdin);
255             inpos = 0;
256             outpos = 0;
257 
258             /* read length table */
259             for(i = 0; i < POLAR_SYMBOLS; i += 2) {
260                 leng_table[i] =     odata[outpos] / 16;
261                 leng_table[i + 1] = odata[outpos] % 16;
262                 outpos++;
263             }
264 
265             /* decode */
266             polar_make_code_table(leng_table, code_table);
267             polar_make_decode_table(leng_table, code_table, decode_table);
268 
269             while(inpos < inlen) {
270                 while(outpos < outlen && code_len < POLAR_MAXLEN) {
271                     code_buf += odata[outpos++] << code_len;
272                     code_len += 8;
273                 }
274                 i = decode_table[code_buf % 65536];
275 
276                 idata[inpos++] = mtf_decode(i);
277                 code_buf >>= leng_table[i];
278                 code_len -=  leng_table[i];
279             }
280             fwrite(idata, 1, inpos, stdout);
281         }
282         return 0;
283     }
284 
285     fprintf(stderr, "%s\n", "Usage: polar [e|d] <input >output");
286     return 0;
287 }

posted on 2012-11-09 19:15  RichSelian  阅读(2842)  评论(0)  编辑  收藏  举报

导航