峰值信噪比PSNR之浮点精度问题

最近程序中需要用到峰值信噪比检验，该算法在计算时需要用到log10运算，经测试发现在Java程序中实现和在Matlab中的实现，针对同一测试数据，可能会有不同的结果。初步推测为不同实现的浮点精度问题。

java代码：

public class Tester {

private static float MSN(Float[] I, Float[] K) {

int m = 0;

float msn = 0;

for (int i = 0; i < I.length; i++) {

if (I[i].isNaN() || K[i].isNaN()) {

} else {

float d = I[i] - K[i];

d *= d;

msn += d;

m++;

}

}

msn = msn / m;

double msnD = Math.sqrt(msn);

msn = new Double(msnD).floatValue();

return msn;

}

/**

* http://en.wikipedia.org/wiki/PSNR
*

* @return Peak signal-to-noise ratio

*/

private static float PSNR(Float[] I, Float[] K) {

float msn = new Double( Math.sqrt(MSN(I, K))).floatValue(), MAX = 255;

float psnr = new Double(20f * Math.log10(MAX / msn)).floatValue();

return psnr;

}

/**

* @param args

*/

public static void main(String[] args) {

// TODO Auto-generated method stub

float[] i = { 2, 3, 4, 5, 6f, 7f, 8f, 9f, 3f, 4f, 6f, 7f, 8f, 10f, 11,

12, 4, 6, 7, 9, 11, 12, 13, 15, 5, 7, 9, 11, 13, 15, 16, 17, 6,

9, 11, 13, 15, 17, 19, 20, 7, 10, 13, 15, 17, 20, 21, 23, 8,

11, 14, 17, 20, 22, 24, 26, 9, 13, 16, 19, 22, 25, 27, 29 };

float[] k = { 3, 3, 4, 6, 7, 8, 9, 10, 3, 4, 5, 6, 8, 9, 10, 11, 4, 5,

6, 8, 10, 12, 13, 14, 6, 7, 8, 10, 12, 15, 16, 17, 7, 8, 10,

13, 15, 18, 20, 21, 9, 10, 12, 15, 18, 21, 23, 24, 10, 11, 14,

17, 20, 23, 25, 26, 11, 12, 14, 17, 21, 24, 26, 28 };

Float[] I = new Float[64];

Float[] K = new Float[64];

int j=0;

for (float ii : i) {

Float iii = new Float(ii);

I[j] = iii;

j++;

}

j=0;

for (float kk : k) {

Float kkk = new Float(kk);

K[j] = kkk;

j++;

}

float psnr = PSNR(I,K);

System.out.println(psnr);

}

}

Matlab代码：

>> A=[2, 3, 4, 5, 6, 7, 8, 9, 3, 4, 6, 7, 8, 10, 11,12, 4, 6, 7, 9, 11, 12, 13, 15, 5, 7, 9, 11, 13, 15, 16, 17, 6,9, 11, 13, 15, 17, 19, 20, 7, 10, 13, 15, 17, 20, 21, 23, 8,11, 14, 17, 20, 22, 24, 26, 9, 13, 16, 19, 22, 25,
27, 29]

A =

Columns 1 through 19

2 3 4 5 6 7 8 9 3 4 6 7 8 10 11 12 4 6 7

Columns 20 through 38

9 11 12 13 15 5 7 9 11 13 15 16 17 6 9 11 13 15 17

Columns 39 through 57

19 20 7 10 13 15 17 20 21 23 8 11 14 17 20 22 24 26 9

Columns 58 through 64

13 16 19 22 25 27 29

>> B=[3, 3, 4, 6, 7, 8, 9, 10, 3, 4, 5, 6, 8, 9, 10, 11, 4, 5,6, 8, 10, 12, 13, 14, 6, 7, 8, 10, 12, 15, 16, 17, 7, 8, 10,13, 15, 18, 20, 21, 9, 10, 12, 15, 18, 21, 23, 24, 10, 11, 14,17, 20, 23, 25, 26, 11, 12, 14, 17, 21, 24, 26, 28]

B =

Columns 1 through 19

3 3 4 6 7 8 9 10 3 4 5 6 8 9 10 11 4 5 6

Columns 20 through 38

8 10 12 13 14 6 7 8 10 12 15 16 17 7 8 10 13 15 18

Columns 39 through 57

20 21 9 10 12 15 18 21 23 24 10 11 14 17 20 23 25 26 11

Columns 58 through 64

12 14 17 21 24 26 28

>> diff = A - B

diff =

Columns 1 through 19

-1 0 0 -1 -1 -1 -1 -1 0 0 1 1 0 1 1 1 0 1 1

Columns 20 through 38

1 1 0 0 1 -1 0 1 1 1 0 0 0 -1 1 1 0 0 -1

Columns 39 through 57

-1 -1 -2 0 1 0 -1 -1 -2 -1 -2 0 0 0 0 -1 -1 0 -2

Columns 58 through 64

1 2 2 1 1 1 1

>> diff_sq = diff.^2

diff_sq =

Columns 1 through 19

1 0 0 1 1 1 1 1 0 0 1 1 0 1 1 1 0 1 1

Columns 20 through 38

1 1 0 0 1 1 0 1 1 1 0 0 0 1 1 1 0 0 1

Columns 39 through 57

1 1 4 0 1 0 1 1 4 1 4 0 0 0 0 1 1 0 4

Columns 58 through 64

1 4 4 1 1 1 1

>> mse = mean(mean(diff_sq));

>> psnr = 10*log10(255^2/mse)

psnr =

48.3393

>> mse = sum(diff.^2)/length(diff)

mse =

0.9531

>> psnr = 10*log10(255^2/mse)

psnr =

48.3393