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