diff options
Diffstat (limited to '2.3-1/src')
-rw-r--r-- | 2.3-1/src/c/linearAlgebra/includes/svd.h | 2 | ||||
-rw-r--r-- | 2.3-1/src/c/linearAlgebra/interfaces/int_sva.h | 2 | ||||
-rw-r--r-- | 2.3-1/src/c/linearAlgebra/interfaces/int_svd.h | 10 | ||||
-rw-r--r-- | 2.3-1/src/c/linearAlgebra/sva/dsvaa.c | 18 | ||||
-rw-r--r-- | 2.3-1/src/c/linearAlgebra/svd/dsvda.c | 82 | ||||
-rw-r--r-- | 2.3-1/src/c/signalProcessing/transforms/dct/cdcta.c | 99 | ||||
-rw-r--r-- | 2.3-1/src/c/signalProcessing/transforms/dct/zdcta.c | 121 | ||||
-rw-r--r-- | 2.3-1/src/c/signalProcessing/transforms/idct/cidcta.c | 45 | ||||
-rw-r--r-- | 2.3-1/src/c/signalProcessing/transforms/idct/zidcta.c | 43 | ||||
-rw-r--r-- | 2.3-1/src/c/string/disp/ddispa.c | 2 |
10 files changed, 317 insertions, 107 deletions
diff --git a/2.3-1/src/c/linearAlgebra/includes/svd.h b/2.3-1/src/c/linearAlgebra/includes/svd.h index 260b87fb..dea681fc 100644 --- a/2.3-1/src/c/linearAlgebra/includes/svd.h +++ b/2.3-1/src/c/linearAlgebra/includes/svd.h @@ -18,7 +18,7 @@ extern "C" { #endif -void dsvda(double *in1,int row,int col,double in2,double nout,double *out1, \ +double dsvda(double tol,double *in1,int row,int col,double in2,double nout,double *out1, \ double *out2,double *out3); void zsvda(doubleComplex *in1,int row,int col,int in2,int nout, doubleComplex *out1,\ doubleComplex *out2,doubleComplex *out3); diff --git a/2.3-1/src/c/linearAlgebra/interfaces/int_sva.h b/2.3-1/src/c/linearAlgebra/interfaces/int_sva.h index 4a2ec56b..f1f8260a 100644 --- a/2.3-1/src/c/linearAlgebra/interfaces/int_sva.h +++ b/2.3-1/src/c/linearAlgebra/interfaces/int_sva.h @@ -18,7 +18,7 @@ extern "C" { #endif #define d2svad2d2d2(in1,size,out1,out2,out3) dsvaa(1,in1,size[0],size[1],0,out1,out2,out3); -#define d2d0svad2d2d2(in1,size1,in2,size2,out1,out2,out3) dsvaa(2,in1,size1[0],size1[1],in2,out1,out2,out3); +#define d2d0svad2d2d2(in1,size1,in2,out1,out2,out3) dsvaa(2,in1,size1[0],size1[1],in2,out1,out2,out3); #ifdef __cplusplus } /* extern "C" */ diff --git a/2.3-1/src/c/linearAlgebra/interfaces/int_svd.h b/2.3-1/src/c/linearAlgebra/interfaces/int_svd.h index 449ee741..8f40bffe 100644 --- a/2.3-1/src/c/linearAlgebra/interfaces/int_svd.h +++ b/2.3-1/src/c/linearAlgebra/interfaces/int_svd.h @@ -17,9 +17,13 @@ extern "C" { #endif -#define d2svdd2(in1,size1,out1) dsvda(in1,size1[0],size1[1],0,1,out1,NULL,NULL) -#define d2g2svdd2d2d2(in1,size1,in2,size2,out1,out2,out3) dsvda(in1,size1[0],size1[1],1,3,out1,out2,out3); -#define d2svdd2d2d2(in1,size1,out1,out2,out3) dsvda(in1,size1[0],size1[1],0,3,out1,out2,out3); +#define d2svdd2(in1,size1,out1) dsvda(0,in1,size1[0],size1[1],0,1,out1,NULL,NULL) +#define d2g2svdd2d2d2(in1,size1,in2,size2,out1,out2,out3) dsvda(0,in1,size1[0],size1[1],1,3,out1,out2,out3); +#define d2svdd2d2d2(in1,size1,out1,out2,out3) dsvda(0,in1,size1[0],size1[1],0,3,out1,out2,out3); + +#define d2svdd2d2d2d0(in1,size1,out1,out2,out3) dsvda(0,in1,size1[0],size1[1],0,4,out1,out2,out3); +#define d2d0svdd2d2d2d0(in1,size1,tol,out1,out2,out3) dsvda(tol,in1,size1[0],size1[1],0,4,out1,out2,out3); + #define z2svdz2(in1,size1,out2) zsvda(in1,size1[0],size1[1],0,1,NULL,out2,NULL); #define z2g2svdz2z2z2(in1,size1,in2,size2,out1,out2,out3) zsvda(in1,size1[0],size1[1],1,3,out1,out2,out3); #define z2svdz2z2z2(in1,size1,out1,out2,out3) zsvda(in1,size1[0],size1[1],0,3,out1,out2,out3); diff --git a/2.3-1/src/c/linearAlgebra/sva/dsvaa.c b/2.3-1/src/c/linearAlgebra/sva/dsvaa.c index ee27eef7..b7d07d8c 100644 --- a/2.3-1/src/c/linearAlgebra/sva/dsvaa.c +++ b/2.3-1/src/c/linearAlgebra/sva/dsvaa.c @@ -33,14 +33,14 @@ void dsvaa(int ninp,double *in1,int row,int col,double in2,double *out1, \ /* Calculation of svd of a given matrix */ double *U,*S,*V; - U = (double *)malloc((double)row*Min(row,col)*sizeof(double)); - S = (double *)malloc((double)Min(row,col)*Min(row,col)*sizeof(double)); - V = (double *)malloc((double)col*Min(row,col)*sizeof(double)); + U = (double *)malloc((double)row*min(row,col)*sizeof(double)); + S = (double *)malloc((double)min(row,col)*min(row,col)*sizeof(double)); + V = (double *)malloc((double)col*min(row,col)*sizeof(double)); - dsvda(in1,M,N,1,3,U,S,V); + dsvda(0,in1,M,N,1,3,U,S,V); if (ninp == 1){ /* [u,s,v] = sva(A) when input is only matrix */ - tol = Max(row,col)*S[0]*eps; + tol = max(row,col)*S[0]*eps; rk = 0; for(i=0;i<col;i++){ if(S[i+i*row] > tol){ @@ -52,7 +52,7 @@ void dsvaa(int ninp,double *in1,int row,int col,double in2,double *out1, \ tol = in2; if(tol > 1){ rk = tol; - if(rk > Min(row,col)){ + if(rk > min(row,col)){ printf("ERROR: Wrong value for input argument !"); out1 = NULL; out2 = NULL; @@ -70,20 +70,20 @@ void dsvaa(int ninp,double *in1,int row,int col,double in2,double *out1, \ } } arow = M; - acol = Min(M,N); + acol = min(M,N); for(i=0;i<arow;i++){ for(j=0;j<rk;j++){ out1[i+j*row]=U[i+j*arow]; } } - arow = Min(M,N); + arow = min(M,N); for(i=0;i<rk;i++){ for(j=0;j<rk;j++){ out2[i+j*(int)rk] = S[i+j*arow]; } } arow = N; - acol = Min(M,N); + acol = min(M,N); for(i=0;i<arow;i++){ for(j=0;j<rk;j++){ out3[i+j*arow] = V[i+j*arow]; diff --git a/2.3-1/src/c/linearAlgebra/svd/dsvda.c b/2.3-1/src/c/linearAlgebra/svd/dsvda.c index e6af3008..c3bcfc29 100644 --- a/2.3-1/src/c/linearAlgebra/svd/dsvda.c +++ b/2.3-1/src/c/linearAlgebra/svd/dsvda.c @@ -27,26 +27,40 @@ int max(int a,int b); extern double dgesvd_(char*,char*,int*,int*,double*,int*,double*,double*,int*,\ double*,int*,double *,int*,int*); +#define eps 2.22044604925e-16 /* pow(2,-52) */ + /* DGESVD computes the singular value decomposition (SVD) of a real M-by-N matrix A, optionally computing the left and/or right singular vectors. The SVD is written A = U * SIGMA * transpose(V) */ -void dsvda(double *in1,int row,int col,double in2,double nout,double *out1, \ - double *out2,double *out3){ +/*Function support - + +s=svd(X) +[U,S,V]=svd(X) +[U,S,V]=svd(X,0) (obsolete) +[U,S,V]=svd(X,"e") +[U,S,V,rk]=svd(X [,tol]) +*/ + +double dsvda(double tol,double *in1,int row,int col,double in2,double nout,double *out1, \ + double *out2,double *out3){ + char JOBU,JOBVT; - int j,k; + int i,j,k; int LDU=1; /*Leading Dimension of U */ int LDVT=1; /*Leading Dimension of VT */ int M = row; int N = col; double *buf; double *S,*U,*VT; - double *WORK; + double *WORK; - if((nout > 1 && in2 == 1) && (M != N)){ /* [U,S,VT] = svd(x,'e') */ + int rk; /*Fourth output if needed */ + + /*if((nout > 1 && in2 == 1) && (M != N)){ // [U,S,VT] = svd(x,'e') if(M > N){ JOBU = 'S'; JOBVT = 'A'; @@ -61,7 +75,7 @@ void dsvda(double *in1,int row,int col,double in2,double nout,double *out1, \ U = (double*) malloc((double) (LDU)*min(M,N)*sizeof(double)); VT = (double*) malloc((double) (LDVT)*N*sizeof(double)); } - else if(nout > 1){ /* [U,S,VT = svd(x)] */ + else */if(nout > 1){ /* [U,S,VT = svd(x)] */ JOBU = 'A'; /*If JOBU = 'A', U contains the M-by-M orthogonal matrix U */ JOBVT = 'A'; /*JOBVT = 'A': all N rows of V**T are returned in the array VT;*/ LDU = M; @@ -74,7 +88,7 @@ void dsvda(double *in1,int row,int col,double in2,double nout,double *out1, \ JOBVT = 'N'; } int LDA = max(1,M); - + /* Making a copy of input matrix */ buf = (double*) malloc((double)M*N*sizeof(double)); memcpy(buf,in1,M*N*sizeof(double)); @@ -100,19 +114,63 @@ void dsvda(double *in1,int row,int col,double in2,double nout,double *out1, \ if(j == k) *((out2+j*(min(M,N)))+k) = *(S+j); else *((out2+j*(min(M,N)))+k) = 0; } - } - dtransposea(VT,LDVT,N,out3); + } + + //dtransposea(VT,LDVT,N,out3); + /*As there is some patch of error in SVD, these lines are added */ + + for(j=1;j<=N;j++){ + for(i=j;i<=N;i++){ + *(out3+i+(j-1)*N-1) = VT[j+(i-1)*N-1]; + *(out3+j+(i-1)*N-1) = VT[i+(j-1)*N-1]; + } + } + /*for(i=0;i<N;i++){ + for(j=0;j<N;j++){ + printf("%lf ",VT[i*row+j]); + } + printf("\n"); + }*/ } else{ - memcpy(out1,U,LDU*min(row,col)*sizeof(double)); + memcpy(out1,U,M*min(M,N)*sizeof(double)); for(j=0;j<min(M,N);j++){ for(k=0;k<min(M,N);k++){ if(j == k) *((out2+j*(min(M,N)))+k) = *(S+j); else *((out2+j*(min(M,N)))+k) = 0; } } - dtransposea(VT,LDVT,N,out3); - } + //dtransposea(VT,LDVT,N,out3); + /*As there is some patch of error in DGESVD, these lines are added */ + /* out3 first taken in some array then will be copied from it. */ + double *outV; + outV = (double *)malloc(N*N*sizeof(double)); + for(j=1;j<=N;j++){ + for(i=j;i<=N;i++){ + *(outV+i+(j-1)*N-1) = VT[j+(i-1)*N-1]; + *(outV+j+(i-1)*N-1) = VT[i+(j-1)*N-1]; + } + } + + for(j=0;j<min(M,N)*N;j++){ + *(out3+j) = *(outV+j); + } + } + + /* From the fortran file of scilab code - if(tol.eq.0.0d0) tol=dble(max(M,N))*eps*stk(lSV) */ + if(tol == 0){ + tol = (double)max(M,N)*eps*S[0]; + } + if(nout == 4){ /*[U,S,VT,rk] = svd(X,tol) where tol - tolerance*/ + rk = 0; + for(i=0;i<min(M,N);i++){ + if(S[i] > tol){ + rk = i+1; + } + } + return rk; + } + return 0; } int min(int a,int b){ diff --git a/2.3-1/src/c/signalProcessing/transforms/dct/cdcta.c b/2.3-1/src/c/signalProcessing/transforms/dct/cdcta.c index 5bc27929..7ff83645 100644 --- a/2.3-1/src/c/signalProcessing/transforms/dct/cdcta.c +++ b/2.3-1/src/c/signalProcessing/transforms/dct/cdcta.c @@ -15,6 +15,7 @@ #include "addition.h" #include "types.h" #include "floatComplex.h" +#include "multiplication.h" /*#include "matrixMultiplication"*/ /*#include <fftw3.h>*/ #include <math.h> @@ -24,10 +25,10 @@ void cdcta(floatComplex *in,int row,int col,int sign,floatComplex *out) int i,j,k,u,v; int n; int x,y; - float res,ress; + float res,ress,vv,ff; float re,z,q,m; - floatComplex accu = DoubleComplex(0, 0); - floatComplex temp,mm; + floatComplex accu = FloatComplex(0, 0); + floatComplex temp,mm,aa,bb,cc; if(sign==-1) { if(row==1) @@ -43,17 +44,25 @@ void cdcta(floatComplex *in,int row,int col,int sign,floatComplex *out) { for(j=0;j<col;j++) { - y=row*j+i; - temp=in[y]*(cos(((M_PI)*(y+1-1./2.)*(x))/n)); + y=row*j+i; + vv = cos(((M_PI)*(y+1-1./2.)*(x))/n); + aa = FloatComplex(vv,0); + temp=cmuls(in[y],aa); out[x]=cadds(out[x],temp); } } if(x==0) - out[x]*=1./(sqrt(n)); + { + vv = 1./(sqrt(n)); + aa = FloatComplex(vv,0); + out[x]=cmuls(out[x],aa); + } else { float res=2./n; - out[x]*=sqrt(res); + res = sqrt(res); + aa = FloatComplex(res,0); + out[x]=cmuls(out[x],aa); } } } @@ -76,29 +85,53 @@ void cdcta(floatComplex *in,int row,int col,int sign,floatComplex *out) y=j*row+i; z=(float)(((float)j+1.0/2.0)*(float)v); q=(float)(M_PI/(float)col); - mm=in[y]*(cos(q*z)); + vv = cos(q*z); + aa = FloatComplex(vv,0); + mm=cmuls(in[y],aa); temp=cadds(temp,mm); } z=(float)(((float)i+1.0/2.0)*(float)u); q=(float)(M_PI/(float)row); - temp*=cos(q*z); + ff = cos(q*z); + bb = FloatComplex(ff,0); + temp=cmuls(temp,bb); out[x]=cadds(out[x],temp); } if(u==0) { - out[x]*=1./sqrt((float)row); + vv = 1./sqrt((float)row); + aa = FloatComplex(vv,0); + out[x]=cmuls(out[x],aa); if(v==0) - out[x]*=1./sqrt((float)col); + { + vv = 1./sqrt((float)col); + aa = FloatComplex(vv,0); + out[x]=cmuls(out[x],aa); + } else - out[x]*=sqrt(2./col); + { + vv = sqrt(2./col); + aa = FloatComplex(vv,0); + out[x]=cmuls(out[x],aa); + } } else { - out[x]*=sqrt(2./row); + vv = sqrt(2./row); + aa = FloatComplex(vv,0); + out[x]=cmuls(out[x],aa); if(v==0) - out[x]*=1./sqrt((float)col); + { + vv = 1./sqrt((float)col); + aa = FloatComplex(vv,0); + out[x]=cmuls(out[x],aa); + } else - out[x]*=sqrt(2./col); + { + vv = sqrt(2./col); + aa = FloatComplex(vv,0); + out[x]=cmuls(out[x],aa); + } } } } @@ -125,12 +158,14 @@ void cdcta(floatComplex *in,int row,int col,int sign,floatComplex *out) if(y==0) { q=res*(cos(((M_PI)*(j)*(v+1./2.))/n)); - out[x]=cadds(out[x],in[y]*q); + aa = FloatComplex(q,0); + out[x]=cadds(out[x],cmuls(in[y],aa)); } else { q=ress*(cos(((M_PI)*(j)*(v+1./2.))/n)); - out[x]=cadds(out[x],in[y]*q); + aa = FloatComplex(q,0); + out[x]=cadds(out[x],cmuls(in[y],aa)); } } } @@ -156,19 +191,37 @@ void cdcta(floatComplex *in,int row,int col,int sign,floatComplex *out) y=row*j+i; mm=in[j*row+i]; z=(float)(((float)v+1.0/2.0)*(float)j); - q=(float)(M_PI/(float)col); - mm=mm*(cos(q*z)); + q=(float)(M_PI/(float)col); + vv = cos(q*z); + aa = FloatComplex(vv,0); + mm=cmuls(mm,aa); if(j==0) - temp=cadds(temp,mm*(1./sqrt((float)col))); + { + vv = 1./sqrt((float)col); + aa = FloatComplex(vv,0); + temp=cadds(temp,cmuls(mm,aa)); + } else - temp=cadds(temp,mm*sqrt(2./col)); + { + vv = sqrt(2./col); + aa = FloatComplex(vv,0); + temp=cadds(temp,cmuls(mm,aa)); + } } z=(float)(((float)u+1.0/2.0)*(float)i); q=(float)(M_PI/(float)row); if(i==0) - out[x]=cadds(out[x],temp*((cos(z*q))*(1./sqrt(row)))); + { + vv = (cos(z*q))*(1./sqrt(row)); + aa = FloatComplex(vv,0); + out[x]=cadds(out[x],cmuls(temp,aa)); + } else - out[x]=cadds(out[x],temp*((cos(z*q))*sqrt(2./row))); + { + vv = (cos(z*q))*sqrt(2./row); + aa = FloatComplex(vv,0); + out[x]=cadds(out[x],cmuls(temp,aa)); + } } } } diff --git a/2.3-1/src/c/signalProcessing/transforms/dct/zdcta.c b/2.3-1/src/c/signalProcessing/transforms/dct/zdcta.c index 0204d682..3ae2e333 100644 --- a/2.3-1/src/c/signalProcessing/transforms/dct/zdcta.c +++ b/2.3-1/src/c/signalProcessing/transforms/dct/zdcta.c @@ -15,6 +15,7 @@ #include "addition.h" #include "types.h" #include "doubleComplex.h" +#include "multiplication.h" /*#include "matrixMultiplication"*/ /*#include <fftw3.h>*/ #include <math.h> @@ -24,10 +25,10 @@ void zdcta(doubleComplex *in,int row,int col,int sign,doubleComplex *out) int i,j,k,u,v; int n; int x,y; - double res,ress; - double re,z,q,m; + double res,ress,vv,ff; + double re,z,q,m; doubleComplex accu = DoubleComplex(0, 0); - doubleComplex temp,mm; + doubleComplex temp,mm,aa,bb,cc; if(sign==-1) { if(row==1) @@ -43,22 +44,30 @@ void zdcta(doubleComplex *in,int row,int col,int sign,doubleComplex *out) { for(j=0;j<col;j++) { - y=row*j+i; - temp=in[y]*(cos(((M_PI)*(y+1-1./2.)*(x))/n)); + y=row*j+i; + vv = cos(((M_PI)*(y+1-1./2.)*(x))/n); + aa = DoubleComplex(vv,0); + temp=zmuls(in[y],aa); out[x]=zadds(out[x],temp); } } if(x==0) - out[x]*=1./(sqrt(n)); + { + vv = 1./(sqrt(n)); + aa = DoubleComplex(vv,0); + out[x]=zmuls(out[x],aa); + } else { - double res=2./n; - out[x]*=sqrt(res); + double res=2./n; + res = sqrt(res); + aa = DoubleComplex(res,0); + out[x]=zmuls(out[x],aa); } } } - } - else + } + else { n=col*row; for(u=0;u<row;u++) @@ -74,31 +83,55 @@ void zdcta(doubleComplex *in,int row,int col,int sign,doubleComplex *out) for(j=0;j<col;j++) { y=j*row+i; - z=(double)(((double)j+1.0/2.0)*(double)v); - q=(double)(M_PI/(double)col); - mm=in[y]*(cos(q*z)); + z=(double )(((double )j+1.0/2.0)*(double )v); + q=(double )(M_PI/(double )col); + vv = cos(q*z); + aa = DoubleComplex(vv,0); + mm=zmuls(in[y],aa); temp=zadds(temp,mm); } - z=(double)(((double)i+1.0/2.0)*(double)u); - q=(double)(M_PI/(double)row); - temp*=cos(q*z); + z=(double )(((double )i+1.0/2.0)*(double )u); + q=(double )(M_PI/(double )row); + ff = cos(q*z); + bb = DoubleComplex(ff,0); + temp=zmuls(temp,bb); out[x]=zadds(out[x],temp); - } + } if(u==0) { - out[x]*=1./sqrt((double)row); + vv = 1./sqrt((double )row); + aa = DoubleComplex(vv,0); + out[x]=zmuls(out[x],aa); if(v==0) - out[x]*=1./sqrt((double)col); + { + vv = 1./sqrt((double )col); + aa = DoubleComplex(vv,0); + out[x]=zmuls(out[x],aa); + } else - out[x]*=sqrt(2./col); + { + vv = sqrt(2./col); + aa = DoubleComplex(vv,0); + out[x]=zmuls(out[x],aa); + } } else { - out[x]*=sqrt(2./row); + vv = sqrt(2./row); + aa = DoubleComplex(vv,0); + out[x]=zmuls(out[x],aa); if(v==0) - out[x]*=1./sqrt((double)col); + { + vv = 1./sqrt((double )col); + aa = DoubleComplex(vv,0); + out[x]=zmuls(out[x],aa); + } else - out[x]*=sqrt(2./col); + { + vv = sqrt(2./col); + aa = DoubleComplex(vv,0); + out[x]=zmuls(out[x],aa); + } } } } @@ -125,12 +158,14 @@ void zdcta(doubleComplex *in,int row,int col,int sign,doubleComplex *out) if(y==0) { q=res*(cos(((M_PI)*(j)*(v+1./2.))/n)); - out[x]=zadds(out[x],in[y]*q); + aa = DoubleComplex(q,0); + out[x]=zadds(out[x],zmuls(in[y],aa)); } else { q=ress*(cos(((M_PI)*(j)*(v+1./2.))/n)); - out[x]=zadds(out[x],in[y]*q); + aa = DoubleComplex(q,0); + out[x]=zadds(out[x],zmuls(in[y],aa)); } } } @@ -155,20 +190,38 @@ void zdcta(doubleComplex *in,int row,int col,int sign,doubleComplex *out) { y=row*j+i; mm=in[j*row+i]; - z=(double)(((double)v+1.0/2.0)*(double)j); - q=(double)(M_PI/(double)col); - mm=mm*(cos(q*z)); + z=(double )(((double )v+1.0/2.0)*(double )j); + q=(double )(M_PI/(double )col); + vv = cos(q*z); + aa = DoubleComplex(vv,0); + mm=zmuls(mm,aa); if(j==0) - temp=zadds(temp,mm*(1./sqrt((double)col))); + { + vv = 1./sqrt((double )col); + aa = DoubleComplex(vv,0); + temp=zadds(temp,zmuls(mm,aa)); + } else - temp=zadds(temp,mm*sqrt(2./col)); + { + vv = sqrt(2./col); + aa = DoubleComplex(vv,0); + temp=zadds(temp,zmuls(mm,aa)); + } } - z=(double)(((double)u+1.0/2.0)*(double)i); - q=(double)(M_PI/(double)row); + z=(double )(((double )u+1.0/2.0)*(double )i); + q=(double )(M_PI/(double )row); if(i==0) - out[x]=zadds(out[x],temp*((cos(z*q))*(1./sqrt(row)))); + { + vv = (cos(z*q))*(1./sqrt(row)); + aa = DoubleComplex(vv,0); + out[x]=zadds(out[x],zmuls(temp,aa)); + } else - out[x]=zadds(out[x],temp*((cos(z*q))*sqrt(2./row))); + { + vv = (cos(z*q))*sqrt(2./row); + aa = DoubleComplex(vv,0); + out[x]=zadds(out[x],zmuls(temp,aa)); + } } } } diff --git a/2.3-1/src/c/signalProcessing/transforms/idct/cidcta.c b/2.3-1/src/c/signalProcessing/transforms/idct/cidcta.c index ec0df0c7..ae98ba19 100644 --- a/2.3-1/src/c/signalProcessing/transforms/idct/cidcta.c +++ b/2.3-1/src/c/signalProcessing/transforms/idct/cidcta.c @@ -14,7 +14,8 @@ #include "idct.h" #include "addition.h" #include "types.h" -#include "floatComplex.h" +#include "floatComplex.h" +#include "multiplication.h" /*#include "matrixMultiplication"*/ /*#include <fftw3.h>*/ #include <math.h> @@ -24,10 +25,10 @@ void cidcta(floatComplex *in,int row,int col,floatComplex *out) int i,j,k,u,v; int n=col; int x,y; - float res,ress; + float res,ress,vv,ff; float re,z,q,m; - floatComplex accu = DoubleComplex(0, 0); - floatComplex temp,mm; + floatComplex accu = FloatComplex(0, 0); + floatComplex temp,mm,aa,bb; if(row==1) { res=1./sqrt(n); @@ -46,12 +47,14 @@ void cidcta(floatComplex *in,int row,int col,floatComplex *out) if(y==0) { q=res*(cos(((M_PI)*(j)*(v+1./2.))/n)); - out[x]=cadds(out[x],in[y]*q); + aa=FloatComplex(q,0); + out[x]=cadds(out[x],cmuls(in[y],aa)); } else { q=ress*(cos(((M_PI)*(j)*(v+1./2.))/n)); - out[x]=cadds(out[x],in[y]*q); + aa=FloatComplex(q,0); + out[x]=cadds(out[x],cmuls(in[y],aa)); } } } @@ -77,19 +80,37 @@ void cidcta(floatComplex *in,int row,int col,floatComplex *out) y=row*j+i; mm=in[j*row+i]; z=(float)(((float)v+1.0/2.0)*(float)j); - q=(float)(M_PI/(float)col); - mm=mm*(cos(q*z)); + q=(float)(M_PI/(float)col); + vv=cos(q*z); + aa=FloatComplex(vv,0); + mm=cmuls(mm,aa); if(j==0) - temp=cadds(temp,mm*(1./sqrt((float)col))); + { + vv=1./sqrt((float)col); + aa=FloatComplex(vv,0); + temp=cadds(temp,cmuls(mm,aa)); + } else - temp=cadds(temp,mm*sqrt(2./col)); + { + vv=sqrt(2./col); + aa=FloatComplex(vv,0); + temp=cadds(temp,cmuls(mm,aa)); + } } z=(float)(((float)u+1.0/2.0)*(float)i); q=(float)(M_PI/(float)row); if(i==0) - out[x]=cadds(out[x],temp*((cos(z*q))*(1./sqrt(row)))); + { + vv=(cos(z*q))*(1./sqrt(row)); + aa=FloatComplex(vv,0); + out[x]=cadds(out[x],cmuls(temp,aa)); + } else - out[x]=cadds(out[x],temp*((cos(z*q))*sqrt(2./row))); + { + vv=(cos(z*q))*sqrt(2./row); + aa=FloatComplex(vv,0); + out[x]=cadds(out[x],cmuls(temp,aa)); + } } } } diff --git a/2.3-1/src/c/signalProcessing/transforms/idct/zidcta.c b/2.3-1/src/c/signalProcessing/transforms/idct/zidcta.c index 2177b18c..cc01c966 100644 --- a/2.3-1/src/c/signalProcessing/transforms/idct/zidcta.c +++ b/2.3-1/src/c/signalProcessing/transforms/idct/zidcta.c @@ -15,6 +15,7 @@ #include "addition.h" #include "types.h" #include "doubleComplex.h" +#include "multiplication.h" /*#include "matrixMultiplication"*/ /*#include <fftw3.h>*/ #include <math.h> @@ -24,10 +25,10 @@ void zidcta(doubleComplex *in,int row,int col,doubleComplex *out) int i,j,k,u,v; int n=col; int x,y; - double res,ress; + double res,ress,vv,ff; double re,z,q,m; doubleComplex accu = DoubleComplex(0, 0); - doubleComplex temp,mm; + doubleComplex temp,mm,aa,bb; if(row==1) { res=1./sqrt(n); @@ -46,17 +47,19 @@ void zidcta(doubleComplex *in,int row,int col,doubleComplex *out) if(y==0) { q=res*(cos(((M_PI)*(j)*(v+1./2.))/n)); - out[x]=zadds(out[x],in[y]*q); + aa=DoubleComplex(q,0); + out[x]=zadds(out[x],zmuls(in[y],aa)); } else { q=ress*(cos(((M_PI)*(j)*(v+1./2.))/n)); - out[x]=zadds(out[x],in[y]*q); + aa=DoubleComplex(q,0); + out[x]=zadds(out[x],zmuls(in[y],aa)); } } } } - + } } else @@ -77,19 +80,37 @@ void zidcta(doubleComplex *in,int row,int col,doubleComplex *out) y=row*j+i; mm=in[j*row+i]; z=(double)(((double)v+1.0/2.0)*(double)j); - q=(double)(M_PI/(double)col); - mm=mm*(cos(q*z)); + q=(double)(M_PI/(double)col); + vv=cos(q*z); + aa=DoubleComplex(vv,0); + mm=zmuls(mm,aa); if(j==0) - temp=zadds(temp,mm*(1./sqrt((double)col))); + { + vv=1./sqrt((double)col); + aa=DoubleComplex(vv,0); + temp=zadds(temp,zmuls(mm,aa)); + } else - temp=zadds(temp,mm*sqrt(2./col)); + { + vv=sqrt(2./col); + aa=DoubleComplex(vv,0); + temp=zadds(temp,zmuls(mm,aa)); + } } z=(double)(((double)u+1.0/2.0)*(double)i); q=(double)(M_PI/(double)row); if(i==0) - out[x]=zadds(out[x],temp*((cos(z*q))*(1./sqrt(row)))); + { + vv=(cos(z*q))*(1./sqrt(row)); + aa=DoubleComplex(vv,0); + out[x]=zadds(out[x],zmuls(temp,aa)); + } else - out[x]=zadds(out[x],temp*((cos(z*q))*sqrt(2./row))); + { + vv=(cos(z*q))*sqrt(2./row); + aa=DoubleComplex(vv,0); + out[x]=zadds(out[x],zmuls(temp,aa)); + } } } } diff --git a/2.3-1/src/c/string/disp/ddispa.c b/2.3-1/src/c/string/disp/ddispa.c index 5e6bb84e..29aaceeb 100644 --- a/2.3-1/src/c/string/disp/ddispa.c +++ b/2.3-1/src/c/string/disp/ddispa.c @@ -16,7 +16,7 @@ double ddispa (double* in, int rows, int columns){ int i = 0,j = 0; for (i = 0; i < rows; ++i) { - for (j=0;j<columns;j++) printf (" %1.20f ", in[i+j*rows]); + for (j=0;j<columns;j++) printf (" %e ", in[i+j*rows]); printf("\n"); } return 0; |