搜索了网上的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”分别处理一下就好。
