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 }