1/* 2 * Copyright (c) 2003 Michael Niedermayer <michaelni@gmx.at> 3 * 4 * This file is part of FFmpeg. 5 * 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. 10 * 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. 15 * 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 19 */ 20 21#include <stdio.h> 22#include <stdlib.h> 23#include <inttypes.h> 24#include <assert.h> 25 26#define FFMIN(a,b) ((a) > (b) ? (b) : (a)) 27#define F 100 28#define SIZE 2048 29 30uint64_t exp16_table[21]={ 31 65537, 32 65538, 33 65540, 34 65544, 35 65552, 36 65568, 37 65600, 38 65664, 39 65793, 40 66050, 41 66568, 42 67616, 43 69763, 44 74262, 45 84150, 46 108051, 47 178145, 48 484249, 49 3578144, 50 195360063, 51 582360139072LL, 52}; 53 54#if 0 55// 16.16 fixpoint exp() 56static unsigned int exp16(unsigned int a){ 57 int i; 58 int out= 1<<16; 59 60 for(i=19;i>=0;i--){ 61 if(a&(1<<i)) 62 out= (out*exp16_table[i] + (1<<15))>>16; 63 } 64 65 return out; 66} 67#endif 68 69// 16.16 fixpoint log() 70static int64_t log16(uint64_t a){ 71 int i; 72 int out=0; 73 74 if(a < 1<<16) 75 return -log16((1LL<<32) / a); 76 a<<=16; 77 78 for(i=20;i>=0;i--){ 79 int64_t b= exp16_table[i]; 80 if(a<(b<<16)) continue; 81 out |= 1<<i; 82 a = ((a/b)<<16) + (((a%b)<<16) + b/2)/b; 83 } 84 return out; 85} 86 87static uint64_t int_sqrt(uint64_t a) 88{ 89 uint64_t ret=0; 90 int s; 91 uint64_t ret_sq=0; 92 93 for(s=31; s>=0; s--){ 94 uint64_t b= ret_sq + (1ULL<<(s*2)) + (ret<<s)*2; 95 if(b<=a){ 96 ret_sq=b; 97 ret+= 1ULL<<s; 98 } 99 } 100 return ret; 101} 102 103int main(int argc,char* argv[]){ 104 int i, j; 105 uint64_t sse=0; 106 uint64_t dev; 107 FILE *f[2]; 108 uint8_t buf[2][SIZE]; 109 uint64_t psnr; 110 int len= argc<4 ? 1 : atoi(argv[3]); 111 int64_t max= (1<<(8*len))-1; 112 int shift= argc<5 ? 0 : atoi(argv[4]); 113 int skip_bytes = argc<6 ? 0 : atoi(argv[5]); 114 int size0=0; 115 int size1=0; 116 117 if(argc<3){ 118 printf("tiny_psnr <file1> <file2> [<elem size> [<shift> [<skip bytes>]]]\n"); 119 printf("For WAV files use the following:\n"); 120 printf("./tiny_psnr file1.wav file2.wav 2 0 44 to skip the header.\n"); 121 return -1; 122 } 123 124 f[0]= fopen(argv[1], "rb"); 125 f[1]= fopen(argv[2], "rb"); 126 if(!f[0] || !f[1]){ 127 fprintf(stderr, "Could not open input files.\n"); 128 return -1; 129 } 130 fseek(f[shift<0], shift < 0 ? -shift : shift, SEEK_SET); 131 132 fseek(f[0],skip_bytes,SEEK_CUR); 133 fseek(f[1],skip_bytes,SEEK_CUR); 134 135 for(;;){ 136 int s0= fread(buf[0], 1, SIZE, f[0]); 137 int s1= fread(buf[1], 1, SIZE, f[1]); 138 139 for(j=0; j<FFMIN(s0,s1); j++){ 140 int64_t a= buf[0][j]; 141 int64_t b= buf[1][j]; 142 if(len==2){ 143 a= (int16_t)(a | (buf[0][++j]<<8)); 144 b= (int16_t)(b | (buf[1][ j]<<8)); 145 } 146 sse += (a-b) * (a-b); 147 } 148 size0 += s0; 149 size1 += s1; 150 if(s0+s1<=0) 151 break; 152 } 153 154 i= FFMIN(size0,size1)/len; 155 if(!i) i=1; 156 dev= int_sqrt( ((sse/i)*F*F) + (((sse%i)*F*F) + i/2)/i ); 157 if(sse) 158 psnr= ((2*log16(max<<16) + log16(i) - log16(sse))*284619LL*F + (1<<31)) / (1LL<<32); 159 else 160 psnr= 1000*F-1; //floating point free infinity :) 161 162 printf("stddev:%5d.%02d PSNR:%3d.%02d bytes:%9d/%9d\n", 163 (int)(dev/F), (int)(dev%F), 164 (int)(psnr/F), (int)(psnr%F), 165 size0, size1); 166 return 0; 167} 168 169 170