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