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 }
原文链接: https://www.cnblogs.com/richselian/archive/2012/11/09/2763162.html
欢迎关注
微信关注下方公众号，第一时间获取干货硬货；公众号内回复【pdf】免费获取数百本计算机经典书籍
原创文章受到原创版权保护。转载请注明出处：https://www.ccppcoding.com/archives/68813
非原创文章文中已经注明原地址，如有侵权，联系删除
关注公众号【高性能架构探索】，第一时间获取最新文章
转载文章受原作者版权保护。转载请注明原作者出处！
