2 * Copyright (c) 2003 Michael Niedermayer <michaelni@gmx.at>
4 * This file is part of FFmpeg.
6 * FFmpeg is free software; you can redistribute it and/or
7 * modify it under the terms of the GNU Lesser General Public
8 * License as published by the Free Software Foundation; either
9 * version 2.1 of the License, or (at your option) any later version.
11 * FFmpeg is distributed in the hope that it will be useful,
12 * but WITHOUT ANY WARRANTY; without even the implied warranty of
13 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
14 * Lesser General Public License for more details.
16 * You should have received a copy of the GNU Lesser General Public
17 * License along with FFmpeg; if not, write to the Free Software
18 * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
29 #include "libavutil/intfloat.h"
30 #include "libavutil/intreadwrite.h"
32 #define FFMIN(a, b) ((a) > (b) ? (b) : (a))
36 uint64_t exp16_table
[21] = {
61 // 16.16 fixpoint exp()
62 static unsigned int exp16(unsigned int a
){
68 out
= (out
*exp16_table
[i
] + (1<<15))>>16;
75 // 16.16 fixpoint log()
76 static int64_t log16(uint64_t a
)
82 return -log16((1LL << 32) / a
);
85 for (i
= 20; i
>= 0; i
--) {
86 int64_t b
= exp16_table
[i
];
90 a
= ((a
/ b
) << 16) + (((a
% b
) << 16) + b
/ 2) / b
;
95 static uint64_t int_sqrt(uint64_t a
)
101 for (s
= 31; s
>= 0; s
--) {
102 uint64_t b
= ret_sq
+ (1ULL << (s
* 2)) + (ret
<< s
) * 2;
111 static int16_t get_s16l(uint8_t *p
)
117 v
.u
= p
[0] | p
[1] << 8;
121 static float get_f32l(uint8_t *p
)
123 union av_intfloat32 v
;
124 v
.i
= p
[0] | p
[1] << 8 | p
[2] << 16 | p
[3] << 24;
128 static double get_f64l(uint8_t *p
)
130 return av_int2double(AV_RL64(p
));
133 static int run_psnr(FILE *f
[2], int len
, int shift
, int skip_bytes
)
138 uint8_t buf
[2][SIZE
];
139 int64_t max
= (1LL << (8 * len
)) - 1;
142 uint64_t maxdist
= 0;
143 double maxdist_d
= 0.0;
146 noseek
= fseek(f
[0], 0, SEEK_SET
) ||
147 fseek(f
[1], 0, SEEK_SET
);
150 for (i
= 0; i
< 2; i
++) {
152 if (fread(p
, 1, 12, f
[i
]) != 12)
154 if (!memcmp(p
, "RIFF", 4) &&
155 !memcmp(p
+ 8, "WAVE", 4)) {
156 if (fread(p
, 1, 8, f
[i
]) != 8)
158 while (memcmp(p
, "data", 4)) {
159 int s
= p
[4] | p
[5] << 8 | p
[6] << 16 | p
[7] << 24;
160 fseek(f
[i
], s
, SEEK_CUR
);
161 if (fread(p
, 1, 8, f
[i
]) != 8)
165 fseek(f
[i
], -12, SEEK_CUR
);
169 fseek(f
[shift
< 0], abs(shift
), SEEK_CUR
);
171 fseek(f
[0], skip_bytes
, SEEK_CUR
);
172 fseek(f
[1], skip_bytes
, SEEK_CUR
);
176 int s0
= fread(buf
[0], 1, SIZE
, f
[0]);
177 int s1
= fread(buf
[1], 1, SIZE
, f
[1]);
179 for (j
= 0; j
< FFMIN(s0
, s1
); j
+= len
) {
186 a
= get_s16l(buf
[0] + j
);
187 b
= get_s16l(buf
[1] + j
);
192 sse
+= (a
- b
) * (a
- b
);
202 a
= get_f64l(buf
[0] + j
);
203 b
= get_f64l(buf
[1] + j
);
205 a
= get_f32l(buf
[0] + j
);
206 b
= get_f32l(buf
[1] + j
);
209 sse_d
+= (a
- b
) * (a
- b
);
210 if (dist
> maxdist_d
)
222 i
= FFMIN(size0
, size1
) / len
;
229 uint64_t dev
= int_sqrt(((sse
/ i
) * F
* F
) + (((sse
% i
) * F
* F
) + i
/ 2) / i
);
231 psnr
= ((2 * log16(max
<< 16) + log16(i
) - log16(sse
)) *
232 284619LL * F
+ (1LL << 31)) / (1LL << 32);
234 psnr
= 1000 * F
- 1; // floating point free infinity :)
236 printf("stddev:%5d.%02d PSNR:%3d.%02d MAXDIFF:%5"PRIu64
" bytes:%9d/%9d\n",
237 (int)(dev
/ F
), (int)(dev
% F
),
238 (int)(psnr
/ F
), (int)(psnr
% F
),
239 maxdist
, size0
, size1
);
245 double psnr
= INT_MAX
;
246 double dev
= sqrt(sse_d
/ i
);
247 uint64_t scale
= (len
== 4) ? (1ULL << 24) : (1ULL << 32);
250 psnr
= 2 * log(DBL_MAX
) - log(i
/ sse_d
);
251 snprintf(psnr_str
, sizeof(psnr_str
), "%5.02f", psnr
);
253 snprintf(psnr_str
, sizeof(psnr_str
), "inf");
255 maxdist
= maxdist_d
* scale
;
257 printf("stddev:%10.2f PSNR:%s MAXDIFF:%10"PRIu64
" bytes:%9d/%9d\n",
258 dev
* scale
, psnr_str
, maxdist
, size0
, size1
);
265 int main(int argc
, char *argv
[])
269 int shift_first
= argc
< 5 ? 0 : atoi(argv
[4]);
270 int skip_bytes
= argc
< 6 ? 0 : atoi(argv
[5]);
271 int shift_last
= shift_first
+ (argc
< 7 ? 0 : atoi(argv
[6]));
274 int max_psnr_shift
= 0;
277 if (!strcmp(argv
[3], "u8")) {
279 } else if (!strcmp(argv
[3], "s16")) {
281 } else if (!strcmp(argv
[3], "f32")) {
283 } else if (!strcmp(argv
[3], "f64")) {
287 len
= strtol(argv
[3], &end
, 0);
288 if (*end
|| len
< 1 || len
> 2) {
289 fprintf(stderr
, "Unsupported sample format: %s\nSupported: u8, s16, f32, f64\n", argv
[3]);
296 printf("tiny_psnr <file1> <file2> [<elem size>|u8|s16|f32|f64 [<shift> [<skip bytes> [<shift search range>]]]]\n");
297 printf("WAV headers are skipped automatically.\n");
301 f
[0] = fopen(argv
[1], "rb");
302 f
[1] = fopen(argv
[2], "rb");
303 if (!f
[0] || !f
[1]) {
304 fprintf(stderr
, "Could not open input files.\n");
308 for (shift
= shift_first
; shift
<= shift_last
; shift
++) {
309 int psnr
= run_psnr(f
, len
, shift
, skip_bytes
);
310 if (psnr
> max_psnr
|| (shift
< 0 && psnr
== max_psnr
)) {
312 max_psnr_shift
= shift
;
315 if (shift_last
> shift_first
)
316 printf("Best PSNR is %3d.%02d for shift %i\n", (int)(max_psnr
/ F
), (int)(max_psnr
% F
), max_psnr_shift
);