搜索了网上的 JPEG 的 matlab 实现方式，发现只有寥寥几个，几乎都是只实现了一半，要么就是哈夫曼编码没有实现，要么就是只算出了哈夫曼的码长计算了一下效率，但是没有实际编码。要么就是太难我看不懂（汗。。。）比如 github 上几位大神的作品。在此只是简单的按照 [1] 中的过程实现了一个简单的 JPEG，没有任何优化，只求简单。在此简述一下 JPEG 编码解码的过程中需要关注的点。
其中仍有许多简化的部分：（1）只计算了一张灰度图的编码解码，如果要是 RGB 通道的图需要额外处理 RGB 到 YUV 的变化，然后再分别对 YUV 进行编码，本文所讨论的即主要对其中的 Y 分量进行了编码。但是这部分的内容网上很容易搜索到。（2）只对固定的图片进行了操作（512*512），因为没有添加 jpeg 文件头，所以很多变量都是我直接指定的。可以添加文件头来增加对多种文件格式的 JPEG。（3）存在一个假设，即每一个 block 在 Z 字形编排后后边全是 0，如果存在一组最后不全是 0，最后一位如果是个数的话这个代码可能就不行了，需要特殊处理一下。(这个可能性很小，但是应该依然存在，我猜。。。)
问题 1：DCT 变换离散余弦变换的实现有三种方式 [1]，第一种是用的矩阵的形式，这个在[4] 中也采用的这种方式。[3]中详细的介绍了 DCT 的原理，非常好评！！！为了方便没有书的同学，在此简述一下方式一：
当进行离散余弦变换时，$ Y = AXA^T $ , 其中 A 即为下边生成的变换矩阵。X 是输入样本矩阵，Y 是变换后的系数矩阵。进行逆变换时：$ X = A^{T}XA $。其中 A 的公式如下：$ A_{ij} = C_i\cos{\frac{(2j+1)i\pi}{2N}} $ 其中：$ C_i = \sqrt{\frac1N} \text{(i=0)}, C_i = \sqrt{\frac2N} \text{(i>0)}$
在 matlab 中只需要一行代码就可以实现：T = dctmtx(8)，当然我们也可以自己创建自己的变换矩阵：
N = 8;
T = zeros(N,N);
for i = 1:N
for j = 1:N
T(i, j) = sqrt(2/N) * cos(((2*j-1)*(i-1)*pi)/(2*N));
end
end
T(1,:) = T(1,:) / sqrt(2);
dctmtx 使用的是矩阵的方法，两种方法结果相似。仅供参考。当然也可以使用 matlab 提供的 dct2()和 idct2()函数。这两个个函数的核心其实是第三种方法。
得到了变换矩阵，进行 DCT 操作。[4]中使用了一种非常简便的方式，BY=blkproc(Y,[8 8],’P1*x*P2′,T,T’); 进行 DCT，可以说非常方便。然后再使用 BY2=blkproc(BY,[8 8],’round(x./P1)’,a); 进行量化。这两个函数使用了 matlab 内置的功能，大大简化了代码。当然出于练习的原因，我还是自己实现了一下这两行代码：
quantization=zeros(X,Y);
for i = 1:8:X
for j = 1:8:Y
mask = input_data(i:i+7,j:j+7);
DCT = T * double(mask) * T’;
quantization(i:i+7,j:j+7) = round(DCT./Luminance_Quantization_Table);
end
end
顺便把量化也做了。在 IDCT 的时候：
data=zeros(X, Y);
for i = 1:8:X
for j = 1:8:Y
mask = decoded_matrix(i:i+7,j:j+7);
mask = mask .* Luminance_Quantization_Table;
data(i:i+7,j:j+7) = T’ * mask * T;
end
end
data = uint8(data);
可以看到恢复的时候有一些损失。注意：一定要把数据转换成 uint8，才能让 imshow()函数生效。
问题 2：Z 字形编排 [4] 中使用了 matlab 自带的函数来处理这个问题，变得异常简单。
% order
order = [1 9 2 3 10 17 25 18 11 4 5 12 19 26 33 …
41 34 27 20 13 6 7 14 21 28 35 42 49 57 50 …
43 36 29 22 15 8 16 23 30 37 44 51 58 59 52 …
45 38 31 24 32 39 46 53 60 61 54 47 40 48 55 …
62 63 56 64];
注意的是 order 的顺序和 Z 字形的序号并不相同，这是为了让这个表来适应 matlab 的特性。使用:
y = im2col(quantization, [8 8], ‘distinct’);
xb = size(y,2);
y = y(order,:);
第一个函数是把整个图划分成 8 *8=64 的若干列，即每一列对应了一个 8 * 8 的块，参数 distinct 是块与块不重叠。若使用默认 sliding 参数即为窗口是滑动的。在 8 * 8 的块排列的时候是按列排列的，窗口移动的过程中也是竖向移动的。（matlab 的默认操作都是以列为基础的，更好的使用列向量的一些特征）
注意：在这里最好把这个生成好的 Y 矩阵保存起来，这样在解码的时候可以和解码后的矩阵进行一下比较，容易发现问题。
恢复的时候，需要使用 order 来生成一个“反 Z 字形”序列：
rev = zeros(1,64);
for k = 1:length(order)
rev(k) = find(order==k);
end
X = 512;
Y = 512;
decoded = decoded(rev,:);
decoded_matrix = col2im(decoded, [8 8], [X Y], ‘distinct’);
问题 3：熵编码（后两个部分几乎都是自己想的，可能效率很低）量化表采用了直接作为变量的形式存了起来。这个在 [3] 中的 github 代码里有 C 语言版的，复制过来稍加改动即可支持 matlab。而对于 Huffman 编码表则是把 [2] 的码表直接以 txt 的形式存了起来。然后使用以下语句:
[ac_RS, ~, ac_code] = textread(‘AC_Hoffman_coding_table.txt’, ‘%s%s%s’);
[dc_RS, ~, dc_code] = textread(‘DC_Hoffman_coding_table.txt’, ‘%s%s%s’);
把 txt 中的内容读进来，读取结果直接是三列，而且 Length 相等，方便下一步操作。当进行转换的时候直接使用类似一种查表的方式（DC 编码方式为例）：
…
if dc==0
% 特殊情况
dc_encode = ’00’;
else
% 先对 SSSS 进行编码
SSSS = floor(log2(abs(dc)))+1;
SSSS_index = strcmp(dc_RS ,string(SSSS));
SSSS_encode = cell2mat(dc_code(SSSS_index));
…
% 对 DIFF 进行编码
if dc > 0
DIFF_encode = dec2bin(dc);
elseif dc < 0
dc_b = abs(dc);
dc_d = bitcmp(uint16(dc_b));
DIFF_encode = dec2bin(dc_d);
DIFF_encode = DIFF_encode(end-SSSS+1:end);
end
注意，在进行小于 0 的数字进行编码时，采取的应该是反码（书上说的是补码，没看懂），这一点在 [3] 中有非常详尽的说明。[3]写的是真的好! 在使用 bitcmp()取反码时，要注意的是取完反码要舍去多余的位，注意 dec2bin 的第二个参数指的是“at least n bits”，是个大坑。所以还是自己截取一下吧。在 AC 编码时几乎和 DC 相同，有几个问题也需要注意，第一个是两种特殊情况的判断。特殊情况 1 时如果后边全是 0 了，就可以直接结束了：
if mask(j:end)==0
ac_encode = [ac_encode ‘1010’];
break
end
特殊情况 2 是超过 15 个 0 了也要特殊处理：
if mask(j)==0
zero_tot = zero_tot + 1;
if zero_tot == 16
ac_encode = [ac_encode ‘11111111001’];
zero_tot = 0;
end
continue
end
然后就是要把 RRRR 和 SSSS 拼接起来，在进行查找：
S = floor(log2(abs(mask(j))))+1;
R_S = [dec2hex(zero_tot) ‘/’ dec2hex(S)];
问题 4：解码解码时先对熵编码进行解码，然后和刚才存好的 y 进行比较，看看是不是一样。然后再恢复 Z 字形，再 IDCT（这两个前边已经写过了）。熵编码解码的时候没次生成一个 64 行一列的 mask，然后拼接到一起。解码的时候先拿着码串取 AC_cod 和 DC_code 中匹配，能匹配成功再进行下一步
index = strcmp(dc_code, tmp);
if any(index)
SSSS = cell2mat(dc_RS(index));
SSSS = str2double(SSSS);
if SSSS == 0
DIFF = ‘0’;
else
DIFF = input_data(p_end+1:p_end+SSSS);
end
使用 any 可以避免一些坑。也要注意 0 的处理。注意首位是 0 的时候，要进行反码操作：
if DIFF(1) == ‘0’
dc_d = bin2dec(DIFF);
dc_b = bitcmp(uint16(dc_d));
dc = -double(mod(dc_b, 2^SSSS));
else
dc = bin2dec(DIFF);
end
dc = -double(mod(dc_b, 2^SSSS)); 这一句如果不加 double 会出现余数为 1 但是取负号就成了 0 的诡异情况。。。难以解释。对 AC 解码时需要注意的是对“0/0”和“F/0”分别处理一下就好。
Reference[1]《多媒体技术基础（第 4 版）》林福宗 清华大学出版社 131-140[2] JPEG 标准推荐的亮度、色度 DC、AC Huffman 编码表[3] JPEG 算法解密[4] 使用 Matlab 实现 JPEG 压缩 (github)[5] Matlab 二进制类型数据相关操作