00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021 #include <stdio.h>
00022 #include <stdlib.h>
00023 #include <inttypes.h>
00024 #include <assert.h>
00025
00026 #define FFMIN(a,b) ((a) > (b) ? (b) : (a))
00027 #define F 100
00028 #define SIZE 2048
00029
00030 uint64_t exp16_table[21]={
00031 65537,
00032 65538,
00033 65540,
00034 65544,
00035 65552,
00036 65568,
00037 65600,
00038 65664,
00039 65793,
00040 66050,
00041 66568,
00042 67616,
00043 69763,
00044 74262,
00045 84150,
00046 108051,
00047 178145,
00048 484249,
00049 3578144,
00050 195360063,
00051 582360139072LL,
00052 };
00053
00054 #if 0
00055
00056 static unsigned int exp16(unsigned int a){
00057 int i;
00058 int out= 1<<16;
00059
00060 for(i=19;i>=0;i--){
00061 if(a&(1<<i))
00062 out= (out*exp16_table[i] + (1<<15))>>16;
00063 }
00064
00065 return out;
00066 }
00067 #endif
00068
00069
00070 static int64_t log16(uint64_t a){
00071 int i;
00072 int out=0;
00073
00074 if(a < 1<<16)
00075 return -log16((1LL<<32) / a);
00076 a<<=16;
00077
00078 for(i=20;i>=0;i--){
00079 int64_t b= exp16_table[i];
00080 if(a<(b<<16)) continue;
00081 out |= 1<<i;
00082 a = ((a/b)<<16) + (((a%b)<<16) + b/2)/b;
00083 }
00084 return out;
00085 }
00086
00087 static uint64_t int_sqrt(uint64_t a)
00088 {
00089 uint64_t ret=0;
00090 int s;
00091 uint64_t ret_sq=0;
00092
00093 for(s=31; s>=0; s--){
00094 uint64_t b= ret_sq + (1ULL<<(s*2)) + (ret<<s)*2;
00095 if(b<=a){
00096 ret_sq=b;
00097 ret+= 1ULL<<s;
00098 }
00099 }
00100 return ret;
00101 }
00102
00103 int main(int argc,char* argv[]){
00104 int i, j;
00105 uint64_t sse=0;
00106 uint64_t dev;
00107 FILE *f[2];
00108 uint8_t buf[2][SIZE];
00109 uint64_t psnr;
00110 int len= argc<4 ? 1 : atoi(argv[3]);
00111 int64_t max= (1<<(8*len))-1;
00112 int shift= argc<5 ? 0 : atoi(argv[4]);
00113 int skip_bytes = argc<6 ? 0 : atoi(argv[5]);
00114 int size0=0;
00115 int size1=0;
00116
00117 if(argc<3){
00118 printf("tiny_psnr <file1> <file2> [<elem size> [<shift> [<skip bytes>]]]\n");
00119 printf("For WAV files use the following:\n");
00120 printf("./tiny_psnr file1.wav file2.wav 2 0 44 to skip the header.\n");
00121 return -1;
00122 }
00123
00124 f[0]= fopen(argv[1], "rb");
00125 f[1]= fopen(argv[2], "rb");
00126 if(!f[0] || !f[1]){
00127 fprintf(stderr, "Could not open input files.\n");
00128 return -1;
00129 }
00130 fseek(f[shift<0], shift < 0 ? -shift : shift, SEEK_SET);
00131
00132 fseek(f[0],skip_bytes,SEEK_CUR);
00133 fseek(f[1],skip_bytes,SEEK_CUR);
00134
00135 for(;;){
00136 int s0= fread(buf[0], 1, SIZE, f[0]);
00137 int s1= fread(buf[1], 1, SIZE, f[1]);
00138
00139 for(j=0; j<FFMIN(s0,s1); j++){
00140 int64_t a= buf[0][j];
00141 int64_t b= buf[1][j];
00142 if(len==2){
00143 a= (int16_t)(a | (buf[0][++j]<<8));
00144 b= (int16_t)(b | (buf[1][ j]<<8));
00145 }
00146 sse += (a-b) * (a-b);
00147 }
00148 size0 += s0;
00149 size1 += s1;
00150 if(s0+s1<=0)
00151 break;
00152 }
00153
00154 i= FFMIN(size0,size1)/len;
00155 if(!i) i=1;
00156 dev= int_sqrt( ((sse/i)*F*F) + (((sse%i)*F*F) + i/2)/i );
00157 if(sse)
00158 psnr= ((2*log16(max<<16) + log16(i) - log16(sse))*284619LL*F + (1<<31)) / (1LL<<32);
00159 else
00160 psnr= 1000*F-1;
00161
00162 printf("stddev:%5d.%02d PSNR:%3d.%02d bytes:%9d/%9d\n",
00163 (int)(dev/F), (int)(dev%F),
00164 (int)(psnr/F), (int)(psnr%F),
00165 size0, size1);
00166 return 0;
00167 }
00168
00169