Mercurial > hg > x
changeset 7:9b1c0825cc77
TFileStream; redefining some function arguments
| author | Wen X <xue.wen@elec.qmul.ac.uk> |
|---|---|
| date | Mon, 11 Oct 2010 17:54:32 +0100 |
| parents | fda5b3561a13 |
| children | f0d2c9b5d3a3 |
| files | hs.cpp quickspec.cpp sinest.cpp sinsyn.cpp sinsyn.h tstream.h |
| diffstat | 6 files changed, 228 insertions(+), 106 deletions(-) [+] |
line wrap: on
line diff
--- a/hs.cpp Wed Oct 06 15:36:50 2010 +0100 +++ b/hs.cpp Mon Oct 11 17:54:32 2010 +0100 @@ -4747,11 +4747,13 @@ CubicSpline(pfr-1, fa, fb, fc, fd, xs, f1, 1, 1); CubicSpline(pfr-1, aa, ab, ac, ad, xs, a1, 1, 1); - for (int fr=0; fr<pfr-1; fr++) Sinusoid(&xrecm[ixs[fr]-dst], 0, offst, aa[fr], ab[fr], ac[fr], ad[fr], fa[fr], fb[fr], fc[fr], fd[fr], p1[fr], p1[fr+1], false); + for (int fr=0; fr<pfr-1; fr++) Sinusoid(offst, &xrecm[ixs[fr]-dst], aa[fr], ab[fr], ac[fr], ad[fr], fa[fr], fb[fr], fc[fr], fd[fr], p1[fr], p1[fr+1], false); // Sinusoid(&xrecm[ixs[0]-dst], -hwid, 0, aa[0], ab[0], ac[0], ad[0], fa[0], fb[0], fc[0], fd[0], p1[0], p1[1], false); // Sinusoid(&xrecm[ixs[pfr-2]-dst], offst, offst+hwid, aa[pfr-2], ab[pfr-2], ac[pfr-2], ad[pfr-2], fa[pfr-2], fb[pfr-2], fc[pfr-2], fd[pfr-2], p1[pfr-2], p1[pfr-1], false); - Sinusoid(&xrecm[ixs[0]-dst], -hwid, 0, 0, 0, 0, ad[0], fa[0], fb[0], fc[0], fd[0], p1[0], p1[1], false); - Sinusoid(&xrecm[ixs[pfr-2]-dst], offst, offst+hwid, 0, 0, 0, ad[pfr-2], fa[pfr-2], fb[pfr-2], fc[pfr-2], fd[pfr-2], p1[pfr-2], p1[pfr-1], false); + double tmpph=p1[0]; Sinusoid(&xrecm[ixs[0]-dst], -hwid, 0, 0, 0, 0, ad[0], fa[0], fb[0], fc[0], fd[0], tmpph, false); + ShiftTrinomial(offst, fa[pfr-1], fb[pfr-1], fc[pfr-1], fd[pfr-1], fa[pfr-2], fb[pfr-2], fc[pfr-2], fd[pfr-2]); + ShiftTrinomial(offst, aa[pfr-1], ab[pfr-1], ac[pfr-1], ad[pfr-1], aa[pfr-2], ab[pfr-2], ac[pfr-2], ad[pfr-2]); + tmpph=p1[pfr-1]; Sinusoid(&xrecm[ixs[pfr-2]-dst], offst, offst+hwid, 0, 0, 0, ad[pfr-1], fa[pfr-1], fb[pfr-1], fc[pfr-1], fd[pfr-1], tmpph, false); if (st_count) {
--- a/quickspec.cpp Wed Oct 06 15:36:50 2010 +0100 +++ b/quickspec.cpp Mon Oct 11 17:54:32 2010 +0100 @@ -1,5 +1,4 @@ //--------------------------------------------------------------------------- - #include <math.h> #include <memory.h> #include "quickspec.h"
--- a/sinest.cpp Wed Oct 06 15:36:50 2010 +0100 +++ b/sinest.cpp Mon Oct 11 17:54:32 2010 +0100 @@ -1625,9 +1625,11 @@ } CubicSpline(Fr-1, fa, fb, fc, fd, xs, ffr, 1, 1); CubicSpline(Fr-1, aa, ab, ac, ad, xs, afr, 1, 1); - for (int fr=0; fr<Fr-1; fr++) Sinusoid(&fs[int(xs[fr])], &as[int(xs[fr])], &phs[int(xs[fr])], &das[int(xs[fr])], 0, Offst, aa[fr], ab[fr], ac[fr], ad[fr], fa[fr], fb[fr], fc[fr], fd[fr], pfr[fr], pfr[fr+1], LogA); - Sinusoid(&fs[int(xs[0])], &as[int(xs[0])], &phs[int(xs[0])], &das[int(xs[0])], -HWid, 0, aa[0], ab[0], ac[0], ad[0], fa[0], fb[0], fc[0], fd[0], pfr[0], pfr[1], LogA); - Sinusoid(&fs[int(xs[Fr-2])], &as[int(xs[Fr-2])], &phs[int(xs[Fr-2])], &das[int(xs[Fr-2])], Offst, Offst+HWid, aa[Fr-2], ab[Fr-2], ac[Fr-2], ad[Fr-2], fa[Fr-2], fb[Fr-2], fc[Fr-2], fd[Fr-2], pfr[Fr-2], pfr[Fr-1], LogA); + for (int fr=0; fr<Fr-1; fr++) Sinusoid(Offst, &fs[int(xs[fr])], &as[int(xs[fr])], &phs[int(xs[fr])], &das[int(xs[fr])], aa[fr], ab[fr], ac[fr], ad[fr], fa[fr], fb[fr], fc[fr], fd[fr], pfr[fr], pfr[fr+1], LogA); + double tmpph=pfr[0]; Sinusoid_direct(&fs[int(xs[0])], &as[int(xs[0])], &phs[int(xs[0])], &das[int(xs[0])], -HWid, 0, aa[0], ab[0], ac[0], ad[0], fa[0], fb[0], fc[0], fd[0], tmpph, LogA); + ShiftTrinomial(xs[Fr-1]-xs[Fr-2], fa[Fr-1], fb[Fr-1], fc[Fr-1], fd[Fr-1], fa[Fr-2], fb[Fr-2], fc[Fr-2], fd[Fr-2]); + ShiftTrinomial(xs[Fr-1]-xs[Fr-2], aa[Fr-1], ab[Fr-1], ac[Fr-1], ad[Fr-1], aa[Fr-2], ab[Fr-2], ac[Fr-2], ad[Fr-2]); + tmpph=pfr[Fr-1]; Sinusoid_direct(&fs[int(xs[Fr-1])], &as[int(xs[Fr-1])], &phs[int(xs[Fr-1])], &das[int(xs[Fr-1])], 0, HWid, aa[Fr-1], ab[Fr-1], ac[Fr-1], ad[Fr-1], fa[Fr-1], fb[Fr-1], fc[Fr-1], fd[Fr-1], tmpph, LogA); delete[] fa; //*/ /* for (int i=0; i<Count; i++) @@ -1716,10 +1718,13 @@ } CubicSpline(Fr-1, fa, fb, fc, fd, xs, ffr, 1, 1); CubicSpline(Fr-1, aa, ab, ac, ad, xs, afr, 1, 1); - for (int fr=0; fr<Fr-1; fr++) Sinusoid(&fs[int(xs[fr])], &as[int(xs[fr])], &phs[int(xs[fr])], &das[int(xs[fr])], 0, Offst, aa[fr], ab[fr], ac[fr], ad[fr], fa[fr], fb[fr], fc[fr], fd[fr], pfr[fr], pfr[fr+1], LogA); - Sinusoid(&fs[int(xs[0])], &as[int(xs[0])], &phs[int(xs[0])], &das[int(xs[0])], -HWid, 0, aa[0], ab[0], ac[0], ad[0], fa[0], fb[0], fc[0], fd[0], pfr[0], pfr[1], LogA); - Sinusoid(&fs[int(xs[Fr-2])], &as[int(xs[Fr-2])], &phs[int(xs[Fr-2])], &das[int(xs[Fr-2])], Offst, Offst+HWid, aa[Fr-2], ab[Fr-2], ac[Fr-2], ad[Fr-2], fa[Fr-2], fb[Fr-2], fc[Fr-2], fd[Fr-2], pfr[Fr-2], pfr[Fr-1], LogA); - delete[] fa; //*/ + for (int fr=0; fr<Fr-1; fr++) Sinusoid(Offst, &fs[int(xs[fr])], &as[int(xs[fr])], &phs[int(xs[fr])], &das[int(xs[fr])], aa[fr], ab[fr], ac[fr], ad[fr], fa[fr], fb[fr], fc[fr], fd[fr], pfr[fr], pfr[fr+1], LogA); + double tmpph=pfr[0]; Sinusoid_direct(&fs[int(xs[0])], &as[int(xs[0])], &phs[int(xs[0])], &das[int(xs[0])], -HWid, 0, aa[0], ab[0], ac[0], ad[0], fa[0], fb[0], fc[0], fd[0], tmpph, LogA); + ShiftTrinomial(xs[Fr-1]-xs[Fr-2], fa[Fr-1], fb[Fr-1], fc[Fr-1], fd[Fr-1], fa[Fr-2], fb[Fr-2], fc[Fr-2], fd[Fr-2]); + ShiftTrinomial(xs[Fr-1]-xs[Fr-2], aa[Fr-1], ab[Fr-1], ac[Fr-1], ad[Fr-1], aa[Fr-2], ab[Fr-2], ac[Fr-2], ad[Fr-2]); + tmpph=pfr[Fr-1]; Sinusoid_direct(&fs[int(xs[Fr-1])], &as[int(xs[Fr-1])], &phs[int(xs[Fr-1])], &das[int(xs[Fr-1])], 0, HWid, aa[Fr-1], ab[Fr-1], ac[Fr-1], ad[Fr-1], fa[Fr-1], fb[Fr-1], fc[Fr-1], fd[Fr-1], tmpph, LogA); + + delete[] fa; //*/ for (int fr=0; fr<Fr; fr++) { @@ -1767,11 +1772,22 @@ int dst=Offst*frst, den=Offst*(fren-1)+Wid, dcount=den-dst; double *f=new double[dcount*2], *ph=&f[dcount]; - for (int fr=frst; fr<fren-1; fr++) Sinusoid(&f[int(xs[fr])-dst], &ph[int(xs[fr])-dst], 0, Offst, fa[fr], fb[fr], fc[fr], fd[fr], phs[fr], phs[fr+1]); - if (frst==0) Sinusoid(&f[int(xs[0])-dst], &ph[int(xs[0])-dst], -HWid, 0, fa[0], fb[0], fc[0], fd[0], phs[0], phs[1]); - else Sinusoid(&f[int(xs[frst-1])-dst], &ph[int(xs[frst-1])-dst], 0, Offst, fa[frst-1], fb[frst-1], fc[frst-1], fd[frst-1], phs[frst-1], phs[frst]); - if (fren==Fr) Sinusoid(&f[int(xs[fren-2])-dst], &ph[int(xs[fren-2])-dst], Offst, Offst+HWid, fa[fren-2], fb[fren-2], fc[fren-2], fd[fren-2], phs[fren-2], phs[fren-1]); - else Sinusoid(&f[int(xs[fren-1])-dst], &ph[int(xs[fren-1])-dst], 0, Offst, fa[fren-1], fb[fren-1], fc[fren-1], fd[fren-1], phs[fren-1], phs[fren]); + for (int fr=frst; fr<fren-1; fr++) Sinusoid(Offst, &f[int(xs[fr])-dst], &ph[int(xs[fr])-dst], fa[fr], fb[fr], fc[fr], fd[fr], phs[fr], phs[fr+1]); + if (frst==0) + { + double tmpph=phs[0]; + Sinusoid_direct(&f[int(xs[0])-dst], &ph[int(xs[0])-dst], -HWid, 0, fa[0], fb[0], fc[0], fd[0], tmpph); + } + else + Sinusoid(Offst, &f[int(xs[frst-1])-dst], &ph[int(xs[frst-1])-dst], fa[frst-1], fb[frst-1], fc[frst-1], fd[frst-1], phs[frst-1], phs[frst]); + if (fren==Fr) + { + double tmpph=phs[Fr-1]; + ShiftTrinomial(Offst, fa[Fr-1], fb[Fr-1], fc[Fr-1], fd[Fr-1], fa[Fr-2], fb[Fr-2], fc[Fr-2], fd[Fr-2]); + Sinusoid_direct(&f[int(xs[Fr-1])-dst], &ph[int(xs[Fr-1])-dst], 0, HWid, fa[Fr-1], fb[Fr-1], fc[Fr-1], fd[Fr-1], tmpph); + } + else + Sinusoid(Offst, &f[int(xs[fren-1])-dst], &ph[int(xs[fren-1])-dst], fa[fren-1], fb[fren-1], fc[fren-1], fd[fren-1], phs[fren-1], phs[fren]); cdouble* y=new cdouble[Wid]; AllocateFFTBuffer(Wid, Amp, W, X); @@ -1924,12 +1940,12 @@ int N=Wid+Offst*(FrCount-1); double* cosine=new double[N], *sine=new double[N]; - Sinusoid(&cosine[hWid], &sine[hWid], -hWid, 0, f3[0], f2[0], f1[0], f0[0], ph); + CosSin(&cosine[hWid], &sine[hWid], -hWid, 0, f3[0], f2[0], f1[0], f0[0], ph); for (int fr=0; fr<FrCount-1; fr++) { int ncentre=hWid+Offst*fr; - if (fr==FrCount-2) Sinusoid(&cosine[ncentre], &sine[ncentre], 0, Wid, f3[fr], f2[fr], f1[fr], f0[fr], ph); - else Sinusoid(&cosine[ncentre], &sine[ncentre], 0, hWid, f3[fr], f2[fr], f1[fr], f0[fr], ph); + if (fr==FrCount-2) CosSin(&cosine[ncentre], &sine[ncentre], 0, Wid, f3[fr], f2[fr], f1[fr], f0[fr], ph); + else CosSin(&cosine[ncentre], &sine[ncentre], 0, hWid, f3[fr], f2[fr], f1[fr], f0[fr], ph); } double err=0; for (int n=0; n<N; n++) {double tmp=cosine[n]-x[n-hWid]; err+=tmp*tmp; tmp=cosine[n]*cosine[n]+sine[n]*sine[n]-1; err+=tmp*tmp;} @@ -1942,13 +1958,13 @@ //store first half of demodulated frame to xs[0:hWid-1] if (fr==0) { - Sinusoid(&refcos[hWid], &refsin[hWid], -hWid, 0, f3[0], f2[0], f1[0], f0[0], ph); + CosSin(&refcos[hWid], &refsin[hWid], -hWid, 0, f3[0], f2[0], f1[0], f0[0], ph); for (int i=0; i<hWid; i++) xs[i].x=ldata[i]*refcos[i], xs[i].y=-ldata[i]*refsin[i]; } else { ph=0; - Sinusoid(refcos, refsin, 0, hWid, f3[fr-1], f2[fr-1], f1[fr-1], f0[fr-1], ph); + CosSin(refcos, refsin, 0, hWid, f3[fr-1], f2[fr-1], f1[fr-1], f0[fr-1], ph); for (int i=0; i<hWid; i++) xs[i].x=ldata[i]*refcos[i], xs[i].y=-ldata[i]*refsin[i]; } @@ -1960,12 +1976,12 @@ //store second half of demodulated frame to xs[hWid:Wid-1] if (fr==FrCount-1) { - Sinusoid(lrefcos, lrefsin, hWid, Wid, f3[FrCount-2], f2[FrCount-2], f1[FrCount-2], f0[FrCount-2], ph); + CosSin(lrefcos, lrefsin, hWid, Wid, f3[FrCount-2], f2[FrCount-2], f1[FrCount-2], f0[FrCount-2], ph); for (int i=hWid; i<Wid; i++) xs[i].x=ldata[i]*lrefcos[i], xs[i].y=-ldata[i]*lrefsin[i]; } else { - Sinusoid(refcos, refsin, 0, hWid, f3[fr], f2[fr], f1[fr], f0[fr], ph); + CosSin(refcos, refsin, 0, hWid, f3[fr], f2[fr], f1[fr], f0[fr], ph); for (int i=hWid; i<Wid; i++) xs[i].x=ldata[i]*lrefcos[i], xs[i].y=-ldata[i]*lrefsin[i]; }
--- a/sinsyn.cpp Wed Oct 06 15:36:50 2010 +0100 +++ b/sinsyn.cpp Mon Oct 11 17:54:32 2010 +0100 @@ -8,6 +8,51 @@ //--------------------------------------------------------------------------- /** + function CosSin: recursive cos-sin generator with trinomial frequency + + In: CountSt, CountEn + f3, f2, f1, f0: trinomial coefficients of frequency + ph: initial phase angle at 0 (NOT at CountSt) + Out: datar[CountSt:CountEn-1], datai[CountSt:CountEn-1]: synthesized pair of cosine and sine functions + ph: phase angle at CountEn + + No return value. +*/ +void CosSin(double* datar, double* datai, int CountSt, int CountEn, double f3, double f2, double f1, double f0, double &ph) +{ + int i; + double dph, ddph, dddph, ddddph, + sph, cph, sdph, cdph, sddph, cddph, sdddph, cdddph, sddddph, cddddph, + p0=ph, p1=2*M_PI*f0, p2=M_PI*f1, p3=2.0*M_PI*f2/3, p4=2.0*M_PI*f3/4, tmp; + if (CountSt==0) + { + dph=p1+p2+p3+p4; ddph=2*p2+6*p3+14*p4; dddph=6*p3+36*p4; ddddph=24*p4; + } + else + { + ph=p0+CountSt*(p1+CountSt*(p2+CountSt*(p3+CountSt*p4))); + dph=p1+p2+p3+p4+CountSt*(2*p2+3*p3+4*p4+CountSt*(3*p3+6*p4+CountSt*4*p4)); + ddph=2*p2+6*p3+14*p4+CountSt*(6*p3+24*p4+CountSt*12*p4); + dddph=6*p3+36*p4+CountSt*24*p4; ddddph=24*p4; + } + sph=sin(ph), cph=cos(ph); + sdph=sin(dph), cdph=cos(dph); + sddph=sin(ddph), cddph=cos(ddph); + sdddph=sin(dddph), cdddph=cos(dddph); + sddddph=sin(ddddph), cddddph=cos(ddddph); + + for (i=CountSt; i<CountEn; i++) + { + datar[i]=cph; datai[i]=sph; + tmp=cph*cdph-sph*sdph; sph=sph*cdph+cph*sdph; cph=tmp; + tmp=cdph*cddph-sdph*sddph; sdph=sdph*cddph+cdph*sddph; cdph=tmp; + tmp=cddph*cdddph-sddph*sdddph; sddph=sddph*cdddph+cddph*sdddph; cddph=tmp; + tmp=cdddph*cddddph-sdddph*sddddph; sdddph=sdddph*cddddph+cdddph*sddddph; cdddph=tmp; + } + ph=p0+CountEn*(p1+CountEn*(p2+CountEn*(p3+CountEn*p4))); +}//CosSin*/ + +/** function Sinuoid: original McAuley-Quatieri synthesizer interpolation between two measurement points. In: T: length from measurement point 1 to measurement point 2 @@ -109,6 +154,70 @@ }//Sinusoid /** + function Sinusoid_direct: synthesizes sinusoid over [CountSt, CountEn) from tronomial coefficients of + amplitude and frequency, direct implementation returning sinusoid coefficients rather than waveform. + + In: CountSt, CountEn + aa, ab, ac, ad: trinomial coefficients of amplitude + fa, fb, fc, fd: trinomial coefficients of frequency + p1: initial phase angle at 0 (NOT at CountSt) + Out: f[CountSt:CountEn-1], a[:], ph[:], da[:]: output buffers. + p1: phase angle at CountEn + + No return value. +*/ +void Sinusoid_direct(double* f, double* a, double* ph, double* da, int CountSt, int CountEn, double aa, double ab, double ac, double ad, + double fa, double fb, double fc, double fd, double &p1, bool LogA) +{ + int i; + if (LogA) + { + for (i=CountSt; i<CountEn; i++) + { + a[i]=exp(ad+i*(ac+i*(ab+i*aa))); + f[i]=fd+i*(fc+i*(fb+i*fa)); + ph[i]=p1+2*M_PI*i*(fd+i*((fc/2)+i*((fb/3)+i*fa/4))); + da[i]=a[i]*(ac+i*(2*ab+3*i*aa)); + } + } + else + { + for (i=CountSt; i<CountEn; i++) + { + a[i]=ad+i*(ac+i*(ab+i*aa)); + f[i]=fd+i*(fc+i*(fb+i*fa)); + ph[i]=p1+2*M_PI*i*(fd+i*((fc/2)+i*((fb/3)+i*fa/4))); + da[i]=ac+i*(2*ab+3*i*aa); + } + } + p1=p1+2*M_PI*i*(fd+i*((fc/2)+i*((fb/3)+i*fa/4))); +}//Sinusoid + +/** + function Sinusoid_direct: synthesizes sinusoid over [CountSt, CountEn) from tronomial coefficients of + frequency, direct implementation returning frequency and phase rather than waveform. + + In: CountSt, CountEn + fa, fb, fc, fd: trinomial coefficients of frequency + p1: initial phase angle at 0 (NOT at CountSt) + Out: f[CountSt:CountEn-1], ph[CountSt:CountEn-1]: output buffers. + p1: phase angle at CountEn + + No return value. +*/ +void Sinusoid_direct(double* f, double* ph, int CountSt, int CountEn, double fa, double fb, double fc, + double fd, double &p1) +{ + int i; + for (i=CountSt; i<CountEn; i++) + { + f[i]=fd+i*(fc+i*(fb+i*fa)); + ph[i]=p1+2*M_PI*i*(fd+i*((fc/2)+i*((fb/3)+i*fa/4))); + } + p1=p1+2*M_PI*i*(fd+i*((fc/2)+i*((fb/3)+i*fa/4))); +}//Sinusoid + +/** function Sinusoid: synthesizes sinusoid over [CountSt, CountEn) from tronomial coefficients of amplitude and frequency, recursive implementation. @@ -294,33 +403,33 @@ function SinusoidExpA: synthesizes complex sinusoid whose log amplitude and frequency are trinomials with phase angle specified at both ends. - In: CountSt, CountEn + In: T a3, a2, a1, a0: trinomial coefficients for log amplitude omg3, omg2, omg1, omg0: trinomial coefficients for angular frequency ph0, ph2: phase angles at 0 and CountEn. add: specifies if the resynthesized sinusoid is to be added to or to replace the content of output buffer - Out: data[CountSt:CountEn-1]: output buffer. + Out: data[T]: output buffer. No return value. */ -void SinusoidExpA(cdouble* data, int CountSt, int CountEn, double a3, double a2, double a1, double a0, +void SinusoidExpA(int T, cdouble* data, double a3, double a2, double a1, double a0, double omg3, double omg2, double omg1, double omg0, double ph0, double ph2, bool add) { double p0=ph0, p1=omg0, p2=0.5*omg1, p3=omg2/3, p4=omg3/4; - double pend=p0+CountEn*(p1+CountEn*(p2+CountEn*(p3+CountEn*p4))); + double pend=p0+T*(p1+T*(p2+T*(p3+T*p4))); int k=floor((pend-ph2)/2/M_PI+0.5); double d=ph2-pend+2*M_PI*k; - double _p=-2*d/CountEn/CountEn/CountEn; - double _q=3*d/CountEn/CountEn; + double _p=-2*d/T/T/T; + double _q=3*d/T/T; - if (add) for (int i=CountSt; i<CountEn; i++) + if (add) for (int i=0; i<T; i++) { double lea=a0+i*(a1+i*(a2+i*a3)); double lph=p0+i*(p1+i*(p2+i*(p3+i*p4))); data[i]+=exp(cdouble(lea, lph+(i*i*(_q+i*_p)))); } - else for (int i=CountSt; i<CountEn; i++) + else for (int i=0; i<T; i++) { double lea=a0+i*(a1+i*(a2+i*a3)); double lph=p0+i*(p1+i*(p2+i*(p3+i*p4))); @@ -474,50 +583,6 @@ ph=p0+CountEn*(p1+CountEn*(p2+CountEn*(p3+CountEn*p4))); } //*/ -/** - function Sinusoid: recursive cos-sin generator with trinomial frequency - - In: CountSt, CountEn - f3, f2, f1, f0: trinomial coefficients of frequency - ph: initial phase angle at 0 (NOT at CountSt) - Out: datar[CountSt:CountEn-1], datai[CountSt:CountEn-1]: synthesized pair of cosine and sine functions - ph: phase angle at CountEn - - No return value. -*/ -void Sinusoid(double* datar, double* datai, int CountSt, int CountEn, double f3, double f2, double f1, double f0, double &ph) -{ - int i; - double dph, ddph, dddph, ddddph, - sph, cph, sdph, cdph, sddph, cddph, sdddph, cdddph, sddddph, cddddph, - p0=ph, p1=2*M_PI*f0, p2=M_PI*f1, p3=2.0*M_PI*f2/3, p4=2.0*M_PI*f3/4, tmp; - if (CountSt==0) - { - dph=p1+p2+p3+p4; ddph=2*p2+6*p3+14*p4; dddph=6*p3+36*p4; ddddph=24*p4; - } - else - { - ph=p0+CountSt*(p1+CountSt*(p2+CountSt*(p3+CountSt*p4))); - dph=p1+p2+p3+p4+CountSt*(2*p2+3*p3+4*p4+CountSt*(3*p3+6*p4+CountSt*4*p4)); - ddph=2*p2+6*p3+14*p4+CountSt*(6*p3+24*p4+CountSt*12*p4); - dddph=6*p3+36*p4+CountSt*24*p4; ddddph=24*p4; - } - sph=sin(ph), cph=cos(ph); - sdph=sin(dph), cdph=cos(dph); - sddph=sin(ddph), cddph=cos(ddph); - sdddph=sin(dddph), cdddph=cos(dddph); - sddddph=sin(ddddph), cddddph=cos(ddddph); - - for (i=CountSt; i<CountEn; i++) - { - datar[i]=cph; datai[i]=sph; - tmp=cph*cdph-sph*sdph; sph=sph*cdph+cph*sdph; cph=tmp; - tmp=cdph*cddph-sdph*sddph; sdph=sdph*cddph+cdph*sddph; cdph=tmp; - tmp=cddph*cdddph-sddph*sdddph; sddph=sddph*cdddph+cddph*sdddph; cddph=tmp; - tmp=cdddph*cddddph-sdddph*sddddph; sdddph=sdddph*cddddph+cdddph*sddddph; cdddph=tmp; - } - ph=p0+CountEn*(p1+CountEn*(p2+CountEn*(p3+CountEn*p4))); -}//Sinusoid*/ /** function Sinusoids: recursive harmonic multi-sinusoid generator @@ -608,24 +673,24 @@ /** function Sinusoid: synthesizes sinusoid piece from trinomial frequency and amplitude coefficients. - In: CountSt, CountEn + In: T aa, ab, ac, ad: trinomial coefficients of amplitude. fa, fb, fc, fd: trinomial coefficients of frequency. ph0, ph2: phase angles at 0 and CountEn. add: specifies if the resynthesized sinusoid is to be added to or to replace the content of output buffer - Out: data[CountSt:CountEn-1]: output buffer. + Out: data[T]: output buffer. No return value. */ -void Sinusoid(double* data, int CountSt, int CountEn, double aa, double ab, double ac, double ad, +void Sinusoid(int T, double* data, double aa, double ab, double ac, double ad, double fa, double fb, double fc, double fd, double ph0, double ph2, bool add) { - double pend=ph0+2*M_PI*CountEn*(fd+CountEn*(fc/2+CountEn*(fb/3+CountEn*fa/4))); + double pend=ph0+2*M_PI*T*(fd+T*(fc/2+T*(fb/3+T*fa/4))); int k=floor((pend-ph2)/2/M_PI+0.5); double d=ph2-pend+2*M_PI*k; - double p=-2*d/CountEn/CountEn/CountEn; - double q=3*d/CountEn/CountEn, a, ph; - for (int i=CountSt; i<CountEn; i++) + double p=-2*d/T/T/T; + double q=3*d/T/T, a, ph; + for (int i=0; i<T; i++) { a=ad+i*(ac+i*(ab+i*aa)); if (a<0) a=0; ph=ph0+2*M_PI*i*(fd+i*((fc/2)+i*((fb/3)+i*fa/4)))+i*i*(q+i*p); @@ -638,7 +703,7 @@ function Sinusoid: synthesizes sinusoid piece from trinomial frequency and amplitude coefficients, returning sinusoid coefficients instead of waveform. - In: CountSt, CountEn + In: T aa, ab, ac, ad: trinomial coefficients of amplitude (or log amplitude if LogA=true) fa, fb, fc, fd: trinomial coefficients of frequency. ph0, ph2: phase angles at 0 and CountEn. @@ -648,22 +713,22 @@ No return value. */ -void Sinusoid(double* f, double* a, double* ph, double* da, int CountSt, int CountEn, double aa, double ab, +void Sinusoid(int T, double* f, double* a, double* ph, double* da, double aa, double ab, double ac, double ad, double fa, double fb, double fc, double fd, double ph0, double ph2, bool LogA) { - double pend=ph0+2*M_PI*CountEn*(fd+CountEn*(fc/2+CountEn*(fb/3+CountEn*fa/4))); + double pend=ph0+2*M_PI*T*(fd+T*(fc/2+T*(fb/3+T*fa/4))); int k=floor((pend-ph2)/2/M_PI+0.5); double d=ph2-pend+2*M_PI*k; - double p=-2*d/CountEn/CountEn/CountEn; - double q=3*d/CountEn/CountEn; - if (LogA) for (int i=CountSt; i<CountEn; i++) + double p=-2*d/T/T/T; + double q=3*d/T/T; + if (LogA) for (int i=0; i<T; i++) { a[i]=exp(ad+i*(ac+i*(ab+i*aa))); if (da) da[i]=a[i]*(ac+i*(2*ab+i*3*aa)); f[i]=fd+i*(fc+i*(fb+i*fa))+i*(2*q+3*i*p)/(2*M_PI); ph[i]=ph0+2*M_PI*i*(fd+i*((fc/2)+i*((fb/3)+i*fa/4)))+i*i*(q+i*p); } - else for (int i=CountSt; i<CountEn; i++) + else for (int i=0; i<T; i++) { a[i]=ad+i*(ac+i*(ab+i*aa)); if (da) da[i]=ac+i*(2*ab+i*3*aa); @@ -675,22 +740,22 @@ /** function Sinusoid: generates trinomial frequency and phase with phase correction. - In: CountSt, CountEn + In: T fa, fb, fc, fd: trinomial coefficients of frequency. ph0, ph2: phase angles at 0 and CountEn. - Out: f[CountSt:CountEn-1], ph[CountSt:CountEn-1]: output buffers holding frequency and phase. + Out: f[T], ph[T]: output buffers holding frequency and phase. No return value. */ -void Sinusoid(double* f, double* ph, int CountSt, int CountEn, double fa, double fb, +void Sinusoid(int T, double* f, double* ph, double fa, double fb, double fc, double fd, double ph0, double ph2) { - double pend=ph0+2*M_PI*CountEn*(fd+CountEn*(fc/2+CountEn*(fb/3+CountEn*fa/4))); + double pend=ph0+2*M_PI*T*(fd+T*(fc/2+T*(fb/3+T*fa/4))); int k=floor((pend-ph2)/2/M_PI+0.5); double d=ph2-pend+2*M_PI*k; - double p=-2*d/CountEn/CountEn/CountEn; - double q=3*d/CountEn/CountEn; - for (int i=CountSt; i<CountEn; i++) + double p=-2*d/T/T/T; + double q=3*d/T/T; + for (int i=0; i<T; i++) { f[i]=fd+i*(fc+i*(fb+i*fa))+i*(2*q+3*i*p)/(2*M_PI); ph[i]=ph0+2*M_PI*i*(fd+i*((fc/2)+i*((fb/3)+i*fa/4)))+i*i*(q+i*p); @@ -765,7 +830,7 @@ *a3=&f3[Fr*4], *a2=&a3[Fr], *a1=&a3[Fr*2], *a0=&a3[Fr*3]; CubicSpline(Fr-1, f3, f2, f1, f0, xs, fs, 1, 1); CubicSpline(Fr-1, a3, a2, a1, a0, xs, as, 1, 1); - for (int fr=0; fr<Fr-1; fr++) Sinusoid(&xrecm[(int)xs[fr]-dst], 0, xs[fr+1]-xs[fr], a3[fr], a2[fr], a1[fr], a0[fr], f3[fr], f2[fr], f1[fr], f0[fr], phs[fr], phs[fr+1], add); + for (int fr=0; fr<Fr-1; fr++) Sinusoid(xs[fr+1]-xs[fr], &xrecm[(int)xs[fr]-dst], a3[fr], a2[fr], a1[fr], a0[fr], f3[fr], f2[fr], f1[fr], f0[fr], phs[fr], phs[fr+1], add); double tmpph=phs[0]; Sinusoid(&xrecm[(int)xs[0]-dst], dst-xs[0], 0, 0, 0, 0, a0[0], f3[0], f2[0], f1[0], f0[0], tmpph, add); //extend the trinomials on [xs[Fr-2], xs[Fr-1]) based at xs[Fr-2] to beyond xs[Fr-1] based at xs[Fr-1]. tmpph=phs[Fr-1];
--- a/sinsyn.h Wed Oct 06 15:36:50 2010 +0100 +++ b/sinsyn.h Mon Oct 11 17:54:32 2010 +0100 @@ -11,13 +11,15 @@ #include "xcomplex.h" //--segmental synthesis routines--------------------------------------------- +void CosSin(double* datar, double* datai, int CountSt, int CountEn, double f3, double f2, double f1, double f0, double &ph); void Sinusoid(double* data, int T, double a1, double a2, double f1, double f2, double p1, double p2, double* a, double* f, double* p, bool ad=true); void Sinusoid(double* data, int T, double a1, double a2, double f1, double f2, double p1, double p2, bool ad=true); -void Sinusoid(double* f, double* a, double* ph, double* da, int CountSt, int CountEn, double aa, double ab, double ac, double ad, double fa, double fb, double fc, double fd, double ph0, double ph2, bool LogA=false); -void Sinusoid(double* f, double* ph, int CountSt, int CountEn, double fa, double fb, double fc, double fd, double ph0, double ph2); -void Sinusoid(double* data, int CountSt, int CountEn, double a3, double a2, double a1, double a0, double f3, double f2, double f1, double f0, double ph0, double ph2, bool add); +void Sinusoid(int T, double* f, double* a, double* ph, double* da, double aa, double ab, double ac, double ad, double fa, double fb, double fc, double fd, double ph0, double ph2, bool LogA=false); +void Sinusoid(int T, double* f, double* ph, double fa, double fb, double fc, double fd, double ph0, double ph2); +void Sinusoid(int T, double* data, double a3, double a2, double a1, double a0, double f3, double f2, double f1, double f0, double ph0, double ph2, bool add); void Sinusoid(double* data, int CountSt, int CountEn, double a3, double a2, double a1, double a0, double f3, double f2, double f1, double f0, double &ph, bool add); -void Sinusoid(double* datar, double* datai, int CountSt, int CountEn, double f3, double f2, double f1, double f0, double &ph); +void Sinusoid_direct(double* f, double* a, double* ph, double* da, int CountSt, int CountEn, double aa, double ab, double ac, double ad, double fa, double fb, double fc, double fd, double &p1, bool LogA=false); +void Sinusoid_direct(double* f, double* ph, int CountSt, int CountEn, double fa, double fb, double fc, double fd, double &p1); void SinusoidExp(cdouble* data, int CountSt, int CountEn, double a3, double a2, double a1, double a0, double omg3, double omg2, double omg1, double omg0, double &ea, double &ph, bool add); void SinusoidExp(int T, cdouble* s, cdouble* ds, int M, cdouble* lamda, double** h, double** dih, cdouble& tmpexp); void SinusoidExp(int T, cdouble* s, int M, cdouble* lamda, double** dih, cdouble& tmpexp); @@ -27,6 +29,7 @@ void SinusoidExpA(int T, cdouble* s, int M, double* p, double* q, double** h, double** dih, double ph1, double ph2); //--multi-segmental synthesis routines--------------------------------------- +void ShiftTrinomial(double T, double& b3, double& b2, double& b1, double& b0, double a3, double a2, double a1, double a0); double* SynthesizeSinusoid(double* xrec, int dst, int den, double* phs, int Fr, double* xs, double* fs, double* as, bool add, bool* terminatetag); double* SynthesizeSinusoidP(double* xrecm, int dst, int den, double* phs, int Fr, double* xs, double* fs, double* as, bool add);
--- a/tstream.h Wed Oct 06 15:36:50 2010 +0100 +++ b/tstream.h Mon Oct 11 17:54:32 2010 +0100 @@ -8,16 +8,53 @@ abstract I/O purposes. */ +#include <stdio.h> + enum TSeekOrigin {soFromBeginning, soFromCurrent, soFromEnd}; class TStream { public: - TStream(); - ~TStream(); - int Read(void*, int){return 0;} - int Write(void*, int){return 0;} - int Seek(int, TSeekOrigin){return Position;} + TStream(){}; + ~TStream(){}; + virtual int Read(void*, int){return 0;} + virtual int Write(void*, int){return 0;} + virtual int Seek(int, TSeekOrigin){return Position;} int Position; }; +enum FileMode {fmRead, fmWrite, fmReadWrite}; +class TFileStream : public TStream +{ +public: + FILE* fp; + TFileStream(char* filename, FileMode mode) : TStream() + { + if (mode==fmWrite) fp=fopen(filename, "wb"); + else if (mode==fmReadWrite) fp=fopen(filename, "rb+"); + else fp=fopen(filename, "rb"); + } + ~TFileStream() + { + fclose(fp); + } + virtual int Read(void* buffer, int size) + { + int result=fread(buffer, 1, size, fp); + Position=ftell(fp); + return result; + } + virtual int Write(void* buffer, int size) + { + int result=fwrite(buffer, 1, size, fp); + Position=ftell(fp); + return result; + } + virtual int Seek(int Offset, TSeekOrigin Origin) + { + fseek(fp, Offset, Origin); + Position=ftell(fp); + return Position; + } +}; + #endif // TSTREAM_H