diff options
Diffstat (limited to 'src/c/signalProcessing')
19 files changed, 715 insertions, 93 deletions
diff --git a/src/c/signalProcessing/%sn/zmodsna.c b/src/c/signalProcessing/%sn/zmodsna.c index 33052a3..85bc0c1 100644 --- a/src/c/signalProcessing/%sn/zmodsna.c +++ b/src/c/signalProcessing/%sn/zmodsna.c @@ -13,12 +13,153 @@ #include<math.h> #include "modsn.h" #include "doubleComplex.h" +#define CA 0.0003 + + +doubleComplex zmodsnsproto(doubleComplex uu,double emmc,doubleComplex* sni) +{ + doubleComplex ans; + double uur,uui; + uur=zreals(uu); + uui=zimags(uu); + double sr,cr,dr; + //Performing Elliptic Function operation for the real values + double a1,b1,c1,d1,emc1,u1; + double em1[14],en1[14]; + int i1,ii1,l1,bo1; + emc1=1-emmc; + u1=uur; + if(emc1) + { + bo1=(emc1<0.0); + if(bo1) + { + d1=1.0-emc1; + emc1/=-1.0/d1; + u1*=(d1=sqrt(d1)); + } + a1=1.0; + dr=1.0; + for(i1=1;i1<=13;i1++) + { + l1=i1; + em1[i1]=a1; + en1[i1]=(emc1=sqrt(emc1)); + c1=0.5*(a1+emc1); + if(fabs(a1-emc1)<=CA*a1)break; + emc1*=a1; + a1=c1; + } + u1*=c1; + sr=sin(u1); + cr=cos(u1); + if(sr) + { + a1=cr/sr; + c1*=a1; + for(ii1=l1;ii1>=1;ii1--) + { + b1=em1[ii1]; + a1*=c1; + c1*=dr; + dr=(en1[ii1]+a1)/(b1+a1); + a1=c1/b1; + } + a1=1.0/sqrt(c1*c1+1.0); + sr=(sr>=0.0?a1:-a1); + cr=c1*(sr); + } + if(bo1) + { + a1=dr; + dr=cr; + cr=a1; + sr/=d1; + } + } + else + { + cr=1.0/cosh(u1); + dr=cr; + sr=tanh(u1); + } + //////////////////////////////////////////////////////////////// + double si,ci,di; + //Performing Elleptic Function operation for the imaginary values + double a,b,c,d,emc,u; + double em[14],en[14]; + int i,ii,l,bo; + //double s1,c1,d1; + emc=emmc; + u=uui; + if(emc) + { + bo=(emc<0.0); + if(bo) + { + d=1.0-emc; + emc/=-1.0/d; + u*=(d=sqrt(d)); + } + a=1.0; + di=1.0; + for(i=1;i<=13;i++) + { + l=i; + em[i]=a; + en[i]=(emc=sqrt(emc)); + c=0.5*(a+emc); + if(fabs(a-emc)<=CA*a)break; + emc*=a; + a=c; + } + u*=c; + si=sin(u); + ci=cos(u); + if(si) + { + a=ci/si; + c*=a; + for(ii=l;ii>=1;ii--) + { + b=em[ii]; + a*=c; + c*=di; + di=(en[ii]+a)/(b+a); + a=c/b; + } + a=1.0/sqrt(c*c+1.0); + si=(si>=0.0?a:-a); + ci=c*(si); + } + if(bo) + { + a=di; + di=ci; + ci=a; + si/=d; + } + } + else + { + ci=1.0/cosh(u); + di=ci; + si=tanh(u); + } + ///////////////////////////////////////////////////////// + double delta; + delta=ci*ci + emmc*sr*sr*si*si; + double snir,snii; + snir=(sr*di)/delta; + snii=(cr*dr*si*ci)/delta; + *sni=DoubleComplex(snir,snii); +} void zmodsna(doubleComplex* uu,int size,double emmc,doubleComplex* sn) { int i; for(i=0;i<size;i++) { - sn[i]=zmodsns(uu[i],emmc); + zmodsnsproto(uu[i],emmc,&sn[i]); } } diff --git a/src/c/signalProcessing/amell/amell.h b/src/c/signalProcessing/amell/amell.h new file mode 100644 index 0000000..30bd6c8 --- /dev/null +++ b/src/c/signalProcessing/amell/amell.h @@ -0,0 +1,27 @@ +/* Copyright (C) 2017 - IIT Bombay - FOSSEE + + This file must be used under the terms of the CeCILL. + This source file is licensed as described in the file COPYING, which + you should have received as part of this distribution. The terms + are also available at + http://www.cecill.info/licences/Licence_CeCILL_V2-en.txt + Author: Ankit Raj + Organization: FOSSEE, IIT Bombay + Email: toolbox@scilab.in + */ + +#ifndef __AMELL_H__ +#define __AMELL_H__ +#include "types.h" + +#ifdef __cplusplus +extern "C" { +#endif + +double damells(double u,double x); + +#ifdef __cplusplus +} /* extern "c" */ +#endif + +#endif /*__AMELL_H__*/ diff --git a/src/c/signalProcessing/amell/damella.c b/src/c/signalProcessing/amell/damella.c new file mode 100644 index 0000000..5c37e2a --- /dev/null +++ b/src/c/signalProcessing/amell/damella.c @@ -0,0 +1,23 @@ +/* Copyright (C) 2017 - IIT Bombay - FOSSEE + + This file must be used under the terms of the CeCILL. + This source file is licensed as described in the file COPYING, which + you should have received as part of this distribution. The terms + are also available at + http://www.cecill.info/licences/Licence_CeCILL_V2-en.txt + Author: Ankit Raj + Organization: FOSSEE, IIT Bombay + Email: toolbox@scilab.in + */ +#include<stdio.h> +#include<math.h> +#include "amell.h" + +void damella(double* u,int size,double x,double* oup) +{ + int i; + for(i=0;i<size;i++) + { + oup[i]=damells(u[i],x); + } +} diff --git a/src/c/signalProcessing/amell/damells.c b/src/c/signalProcessing/amell/damells.c new file mode 100644 index 0000000..90c2053 --- /dev/null +++ b/src/c/signalProcessing/amell/damells.c @@ -0,0 +1,57 @@ +/* Copyright (C) 2017 - IIT Bombay - FOSSEE + + This file must be used under the terms of the CeCILL. + This source file is licensed as described in the file COPYING, which + you should have received as part of this distribution. The terms + are also available at + http://www.cecill.info/licences/Licence_CeCILL_V2-en.txt + Author: Ankit Raj + Organization: FOSSEE, IIT Bombay + Email: toolbox@scilab.in + */ + +#include<stdio.h> +#include<math.h> +#include "amell.h" +#define N 30 +#define DBL_EPSILON 2.2204460492503131E-16 + + +double damells(double u,double x) +{ + double a[N+1]; + double g[N+1]; + double c[N+1]; + double two_n; + double phi; + double k; + int n; + k=(long double)fabs(x); + if(k==1.0) + return 0; + if(k>1.0) + printf("Wrong type of input argument type #2"); + + a[0]=1.0; + g[0]=sqrt(1.0-k*k); + c[0]=k; + two_n=1.0; + for(n=0;n<N;n++) + { + if(fabs(a[n]-g[n])<(a[n]*DBL_EPSILON)) + { + break; + } + two_n+=two_n; + a[n+1]=0.5*(a[n]+g[n]); + g[n+1]=sqrt(a[n]*g[n]); + c[n+1]=0.5*(a[n]-g[n]); + } + phi=two_n*a[n]*u; + for(;n>0;n--) + { + phi=0.5*(phi+asin(c[n]*sin(phi)/a[n])); + } + return (double)phi; +} + diff --git a/src/c/signalProcessing/amell/int_amell.h b/src/c/signalProcessing/amell/int_amell.h new file mode 100644 index 0000000..5d0c86f --- /dev/null +++ b/src/c/signalProcessing/amell/int_amell.h @@ -0,0 +1,18 @@ +/* Copyright (C) 2017 - IIT Bombay - FOSSEE + + This file must be used under the terms of the CeCILL. + This source file is licensed as described in the file COPYING, which + you should have received as part of this distribution. The terms + are also available at + http://www.cecill.info/licences/Licence_CeCILL_V2-en.txt + Author: Ankit Raj + Organization: FOSSEE, IIT Bombay + Email: toolbox@scilab.in + */ + +#ifndef __INT_AMELL_H__ +#define __INT_AMELL_H__ + +#define d0d0amelld0(u,x) damells(u,x) + +#endif /* !__INT_AMELL_H__! */ diff --git a/src/c/signalProcessing/ell1mag/dell1maga.c b/src/c/signalProcessing/ell1mag/dell1maga.c new file mode 100644 index 0000000..9af0c8e --- /dev/null +++ b/src/c/signalProcessing/ell1mag/dell1maga.c @@ -0,0 +1,40 @@ +/* Copyright (C) 2017 - IIT Bombay - FOSSEE + + This file must be used under the terms of the CeCILL. + This source file is licensed as described in the file COPYING, which + you should have received as part of this distribution. The terms + are also available at + http://www.cecill.info/licences/Licence_CeCILL_V2-en.txt + Author: Ankit Raj + Organization: FOSSEE, IIT Bombay + Email: toolbox@scilab.in + */ + +#include<stdio.h> +#include "modsn.h" +#include "ell1mag.h" +#include "doubleComplex.h" + +void dell1maga(double eps,double m1,double* z,int size,double* oup ) +{ + double s[size]; + int i; + for(i=0;i<size;i++) + { + s[i]=zmodsns(z[i],m1); + } + double un[size]; + int j; + for(j=0;j<size;j++) + { + un[j]=1; + } + double v; + int k; + for(k=0;k<size;k++) + { + v=un[k]/(un[k]+(eps*eps*s[k]*s[k])); + oup[k]=v; + } +} + diff --git a/src/c/signalProcessing/ell1mag/ell1mag.h b/src/c/signalProcessing/ell1mag/ell1mag.h new file mode 100644 index 0000000..8fffc0c --- /dev/null +++ b/src/c/signalProcessing/ell1mag/ell1mag.h @@ -0,0 +1,28 @@ +/* Copyright (C) 2017 - IIT Bombay - FOSSEE + + This file must be used under the terms of the CeCILL. + This source file is licensed as described in the file COPYING, which + you should have received as part of this distribution. The terms + are also available at + http://www.cecill.info/licences/Licence_CeCILL_V2-en.txt + Author: Ankit Raj + Organization: FOSSEE, IIT Bombay + Email: toolbox@scilab.in + */ + +#ifndef __ELL1MAG_H__ +#define __ELL1MAG_H__ +#include "types.h" +#include "doubleComplex.h" + +#ifdef __cplusplus +extern "C" { +#endif + +void dell1mags(double eps,double m1,doubleComplex* z,int size,double* oup); + +#ifdef __cplusplus +} /* extern "C" */ +#endif + +#endif /*__ELL1MAG_H__*/ diff --git a/src/c/signalProcessing/ell1mag/int_ell1mag.h b/src/c/signalProcessing/ell1mag/int_ell1mag.h new file mode 100644 index 0000000..590a0ab --- /dev/null +++ b/src/c/signalProcessing/ell1mag/int_ell1mag.h @@ -0,0 +1,18 @@ +/* Copyright (C) 2017 - IIT Bombay - FOSSEE + + This file must be used under the terms of the CeCILL. + This source file is licensed as described in the file COPYING, which + you should have received as part of this distribution. The terms + are also available at + http://www.cecill.info/licences/Licence_CeCILL_V2-en.txt + Author: Ankit Raj + Organization: FOSSEE, IIT Bombay + Email: toolbox@scilab.in + */ + +#ifndef __INT_ELL1MAG_H__ +#define __INT_ELL1MAG_H__ + +#define d0d0z2ell1magd2(eps,m1,z,size,oup) dell1mags(eps,m1,z,size,oup) + +#endif /* !__INT_ELL1MAG_H__! */ diff --git a/src/c/signalProcessing/ell1mag/zell1maga.c b/src/c/signalProcessing/ell1mag/zell1maga.c new file mode 100644 index 0000000..6e7a6f9 --- /dev/null +++ b/src/c/signalProcessing/ell1mag/zell1maga.c @@ -0,0 +1,40 @@ +/* Copyright (C) 2017 - IIT Bombay - FOSSEE + + This file must be used under the terms of the CeCILL. + This source file is licensed as described in the file COPYING, which + you should have received as part of this distribution. The terms + are also available at + http://www.cecill.info/licences/Licence_CeCILL_V2-en.txt + Author: Ankit Raj + Organization: FOSSEE, IIT Bombay + Email: toolbox@scilab.in + */ + +#include<stdio.h> +#include "modsn.h" +#include "ell1mag.h" +#include "doubleComplex.h" + +void zell1maga(double eps,double m1,doubleComplex* z,int size,double* oup ) +{ + doubleComplex s[size]; + int i; + for(i=0;i<size;i++) + { + s[i]=zmodsns(z[i],m1); + } + double un[size]; + int j; + for(j=0;j<size;j++) + { + un[j]=1; + } + doubleComplex v; + int k; + for(k=0;k<size;k++) + { + v=un[k]/(un[k]+(eps*eps*s[k]*s[k])); + oup[k]=zreals(v); + } +} + diff --git a/src/c/signalProcessing/filt_sinc/dfilt_sincs.c b/src/c/signalProcessing/filt_sinc/dfilt_sincs.c index 65aaaa6..1b7d1b1 100644 --- a/src/c/signalProcessing/filt_sinc/dfilt_sincs.c +++ b/src/c/signalProcessing/filt_sinc/dfilt_sincs.c @@ -46,13 +46,4 @@ void dfilt_sincs(double N,double fc,double* oup) oup[k]=xn[k]/xd[k]; } } -/* -int main() -{ - int n; - double fl; - n=5; - fl=0.2; - filt_sinc(n,fl); -} -*/ + diff --git a/src/c/signalProcessing/includes/amell.h b/src/c/signalProcessing/includes/amell.h new file mode 100644 index 0000000..2336d3c --- /dev/null +++ b/src/c/signalProcessing/includes/amell.h @@ -0,0 +1,28 @@ +/* Copyright (C) 2017 - IIT Bombay - FOSSEE + + This file must be used under the terms of the CeCILL. + This source file is licensed as described in the file COPYING, which + you should have received as part of this distribution. The terms + are also available at + http://www.cecill.info/licences/Licence_CeCILL_V2-en.txt + Author: Ankit Raj + Organization: FOSSEE, IIT Bombay + Email: toolbox@scilab.in + */ + +#ifndef __AMELL_H__ +#define __AMELL_H__ +#include "types.h" + +#ifdef __cplusplus +extern "C" { +#endif + +double damells(double u,double x); +void damella(double* u,int size,double x,double* oup); + +#ifdef __cplusplus +} /* extern "c" */ +#endif + +#endif /*__AMELL_H__*/ diff --git a/src/c/signalProcessing/includes/ell1mag.h b/src/c/signalProcessing/includes/ell1mag.h new file mode 100644 index 0000000..e881cca --- /dev/null +++ b/src/c/signalProcessing/includes/ell1mag.h @@ -0,0 +1,29 @@ +/* Copyright (C) 2017 - IIT Bombay - FOSSEE + + This file must be used under the terms of the CeCILL. + This source file is licensed as described in the file COPYING, which + you should have received as part of this distribution. The terms + are also available at + http://www.cecill.info/licences/Licence_CeCILL_V2-en.txt + Author: Ankit Raj + Organization: FOSSEE, IIT Bombay + Email: toolbox@scilab.in + */ + +#ifndef __ELL1MAG_H__ +#define __ELL1MAG_H__ +#include "types.h" +#include "doubleComplex.h" + +#ifdef __cplusplus +extern "C" { +#endif + +void zell1maga(double eps,double m1,doubleComplex* z,int size,double* oup); +void dell1maga(double eps,double m1,double* z,int size,double* oup); + +#ifdef __cplusplus +} /* extern "C" */ +#endif + +#endif /*__ELL1MAG_H__*/ diff --git a/src/c/signalProcessing/interfaces/int_amell.h b/src/c/signalProcessing/interfaces/int_amell.h new file mode 100644 index 0000000..10719ac --- /dev/null +++ b/src/c/signalProcessing/interfaces/int_amell.h @@ -0,0 +1,19 @@ +/* Copyright (C) 2017 - IIT Bombay - FOSSEE + + This file must be used under the terms of the CeCILL. + This source file is licensed as described in the file COPYING, which + you should have received as part of this distribution. The terms + are also available at + http://www.cecill.info/licences/Licence_CeCILL_V2-en.txt + Author: Ankit Raj + Organization: FOSSEE, IIT Bombay + Email: toolbox@scilab.in + */ + +#ifndef __INT_AMELL_H__ +#define __INT_AMELL_H__ + +#define d0d0amelld0(u,x) damells(u,x) +#define d2d0amelld2(u,size,x,oup) damella(u,size[1],x,oup) + +#endif /* !__INT_AMELL_H__! */ diff --git a/src/c/signalProcessing/interfaces/int_ell1mag.h b/src/c/signalProcessing/interfaces/int_ell1mag.h new file mode 100644 index 0000000..c30ffef --- /dev/null +++ b/src/c/signalProcessing/interfaces/int_ell1mag.h @@ -0,0 +1,19 @@ +/* Copyright (C) 2017 - IIT Bombay - FOSSEE + + This file must be used under the terms of the CeCILL. + This source file is licensed as described in the file COPYING, which + you should have received as part of this distribution. The terms + are also available at + http://www.cecill.info/licences/Licence_CeCILL_V2-en.txt + Author: Ankit Raj + Organization: FOSSEE, IIT Bombay + Email: toolbox@scilab.in + */ + +#ifndef __INT_ELL1MAG_H__ +#define __INT_ELL1MAG_H__ + +#define d0d0z2ell1magd2(eps,m1,z,size,oup) zell1maga(eps,m1,z,size[1],oup) +#define d0d0d2ell1magd2(eps,m1,z,size,oup) dell1maga(eps,m1,z,size[1],oup) + +#endif /* !__INT_ELL1MAG_H__! */ diff --git a/src/c/signalProcessing/interfaces/int_modsn.h b/src/c/signalProcessing/interfaces/int_modsn.h index 0d32eb0..56c8f8c 100644 --- a/src/c/signalProcessing/interfaces/int_modsn.h +++ b/src/c/signalProcessing/interfaces/int_modsn.h @@ -18,4 +18,4 @@ #define d2d0modsnd2(uu,size,emmc,sn) dmodsna(uu,size[1],emmc,sn) #define z2d0modsnz2(uu,size,emmc,sn) zmodsna(uu,size[1],emmc,sn) -#endif /* !INT_MODSN_H__! */ +#endif /* !__INT_MODSN_H__! */ diff --git a/src/c/signalProcessing/transforms/dct/cdcta.c b/src/c/signalProcessing/transforms/dct/cdcta.c index 7a006a0..7ff8364 100644 --- a/src/c/signalProcessing/transforms/dct/cdcta.c +++ b/src/c/signalProcessing/transforms/dct/cdcta.c @@ -25,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 = FloatComplex(0, 0); - floatComplex temp,mm; + floatComplex temp,mm,aa,bb,cc; if(sign==-1) { if(row==1) @@ -44,20 +44,25 @@ void cdcta(floatComplex *in,int row,int col,int sign,floatComplex *out) { for(j=0;j<col;j++) { - y=row*j+i; - float a; - a=(cos(((M_PI)*(y+1-1./2.)*(x))/n)); - floatComplex b=FloatComplex(a,0); - temp=cmuls(in[y],b); + 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); } } } @@ -80,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); + } } } } @@ -129,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)); } } } @@ -160,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/src/c/signalProcessing/transforms/dct/zdcta.c b/src/c/signalProcessing/transforms/dct/zdcta.c index 0204d68..3ae2e33 100644 --- a/src/c/signalProcessing/transforms/dct/zdcta.c +++ b/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/src/c/signalProcessing/transforms/idct/cidcta.c b/src/c/signalProcessing/transforms/idct/cidcta.c index e6c746c..ae98ba1 100644 --- a/src/c/signalProcessing/transforms/idct/cidcta.c +++ b/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 = FloatComplex(0, 0); - floatComplex temp,mm; + 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/src/c/signalProcessing/transforms/idct/zidcta.c b/src/c/signalProcessing/transforms/idct/zidcta.c index 2177b18..cc01c96 100644 --- a/src/c/signalProcessing/transforms/idct/zidcta.c +++ b/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)); + } } } } |