ref: 781b09f7c290edce02c9169d1c67986e46ad4f2f
parent: 921fc834fa7d87e16b1be27fc7f808fc1ad50a69
author: menno <menno>
date: Thu Dec 23 09:34:51 EST 1999
New fft in psycho model (fast) and new unpredictabilty measure calculation
--- a/faac.dsp
+++ b/faac.dsp
@@ -51,6 +51,7 @@
LINK32=link.exe
# ADD BASE LINK32 kernel32.lib user32.lib gdi32.lib winspool.lib comdlg32.lib advapi32.lib shell32.lib ole32.lib oleaut32.lib uuid.lib odbc32.lib odbccp32.lib kernel32.lib user32.lib gdi32.lib winspool.lib comdlg32.lib advapi32.lib shell32.lib ole32.lib oleaut32.lib uuid.lib odbc32.lib odbccp32.lib /nologo /subsystem:console /machine:I386
# ADD LINK32 libsndfile.lib kernel32.lib user32.lib gdi32.lib winspool.lib comdlg32.lib advapi32.lib shell32.lib ole32.lib oleaut32.lib uuid.lib odbc32.lib odbccp32.lib /nologo /subsystem:console /machine:I386
+# SUBTRACT LINK32 /profile
!ELSEIF "$(CFG)" == "faac - Win32 Debug"
@@ -74,7 +75,7 @@
# ADD BSC32 /nologo
LINK32=link.exe
# ADD BASE LINK32 kernel32.lib user32.lib gdi32.lib winspool.lib comdlg32.lib advapi32.lib shell32.lib ole32.lib oleaut32.lib uuid.lib odbc32.lib odbccp32.lib kernel32.lib user32.lib gdi32.lib winspool.lib comdlg32.lib advapi32.lib shell32.lib ole32.lib oleaut32.lib uuid.lib odbc32.lib odbccp32.lib /nologo /subsystem:console /debug /machine:I386 /pdbtype:sept
-# ADD LINK32 libsndfile.lib kernel32.lib user32.lib gdi32.lib winspool.lib comdlg32.lib advapi32.lib shell32.lib ole32.lib oleaut32.lib uuid.lib odbc32.lib odbccp32.lib /nologo /subsystem:console /debug /machine:I386 /pdbtype:sept
+# ADD LINK32 libsndfile.lib kernel32.lib user32.lib gdi32.lib winspool.lib comdlg32.lib advapi32.lib shell32.lib ole32.lib oleaut32.lib uuid.lib odbc32.lib odbccp32.lib /nologo /subsystem:console /profile /debug /machine:I386
!ENDIF
@@ -111,6 +112,14 @@
# End Source File
# Begin Source File
+SOURCE=.\fastfft.c
+# End Source File
+# Begin Source File
+
+SOURCE=.\fft.c
+# End Source File
+# Begin Source File
+
SOURCE=.\imdct.c
# End Source File
# Begin Source File
@@ -140,10 +149,6 @@
# Begin Source File
SOURCE=.\transfo.c
-# End Source File
-# Begin Source File
-
-SOURCE=.\fastfft.c
# End Source File
# End Group
# End Target
--- a/faac_dll.dsp
+++ b/faac_dll.dsp
@@ -43,7 +43,7 @@
# PROP Ignore_Export_Lib 0
# PROP Target_Dir ""
# ADD BASE CPP /nologo /MT /W3 /GX /O2 /D "WIN32" /D "NDEBUG" /D "_WINDOWS" /D "_MBCS" /D "_USRDLL" /D "FAAC_DLL_EXPORTS" /YX /FD /c
-# ADD CPP /nologo /MT /W1 /O2 /D "NDEBUG" /D "FAAC_DLL" /D "WIN32" /D "_WINDOWS" /D "_MBCS" /D "_USRDLL" /D "FAAC_DLL_EXPORTS" /YX /FD /c
+# ADD CPP /nologo /MT /O2 /D "NDEBUG" /D "FAAC_DLL" /D "WIN32" /D "_WINDOWS" /D "_MBCS" /D "_USRDLL" /D "FAAC_DLL_EXPORTS" /YX /FD /c
# ADD BASE MTL /nologo /D "NDEBUG" /mktyplib203 /win32
# ADD MTL /nologo /D "NDEBUG" /mktyplib203 /win32
# ADD BASE RSC /l 0x413 /d "NDEBUG"
@@ -116,6 +116,14 @@
# End Source File
# Begin Source File
+SOURCE=.\fastfft.c
+# End Source File
+# Begin Source File
+
+SOURCE=.\fft.c
+# End Source File
+# Begin Source File
+
SOURCE=.\imdct.c
# End Source File
# Begin Source File
@@ -145,10 +153,6 @@
# Begin Source File
SOURCE=.\transfo.c
-# End Source File
-# Begin Source File
-
-SOURCE=.\fastfft.c
# End Source File
# End Group
# Begin Group "Resource Files"
--- /dev/null
+++ b/fft.c
@@ -1,0 +1,203 @@
+/*
+** FFT and FHT routines
+** Copyright 1988, 1993; Ron Mayer
+**
+** fht(fz,n);
+** Does a hartley transform of "n" points in the array "fz".
+**
+** NOTE: This routine uses at least 2 patented algorithms, and may be
+** under the restrictions of a bunch of different organizations.
+** Although I wrote it completely myself; it is kind of a derivative
+** of a routine I once authored and released under the GPL, so it
+** may fall under the free software foundation's restrictions;
+** it was worked on as a Stanford Univ project, so they claim
+** some rights to it; it was further optimized at work here, so
+** I think this company claims parts of it. The patents are
+** held by R. Bracewell (the FHT algorithm) and O. Buneman (the
+** trig generator), both at Stanford Univ.
+** If it were up to me, I'd say go do whatever you want with it;
+** but it would be polite to give credit to the following people
+** if you use this anywhere:
+** Euler - probable inventor of the fourier transform.
+** Gauss - probable inventor of the FFT.
+** Hartley - probable inventor of the hartley transform.
+** Buneman - for a really cool trig generator
+** Mayer(me) - for authoring this particular version and
+** including all the optimizations in one package.
+** Thanks,
+** Ron Mayer; mayer@acuson.com
+** and added some optimization by
+** Mather - idea of using lookup table
+** Takehiro - some dirty hack for speed up
+*/
+
+#include <math.h>
+
+#define SQRT2 sqrt(2)
+
+static float costab[20]=
+ {+ .00000000000000000000000000000000000000000000000000,
+ .70710678118654752440084436210484903928483593768847,
+ .92387953251128675612818318939678828682241662586364,
+ .98078528040323044912618223613423903697393373089333,
+ .99518472667219688624483695310947992157547486872985,
+ .99879545620517239271477160475910069444320361470461,
+ .99969881869620422011576564966617219685006108125772,
+ .99992470183914454092164649119638322435060646880221,
+ .99998117528260114265699043772856771617391725094433,
+ .99999529380957617151158012570011989955298763362218,
+ .99999882345170190992902571017152601904826792288976,
+ .99999970586288221916022821773876567711626389934930,
+ .99999992646571785114473148070738785694820115568892,
+ .99999998161642929380834691540290971450507605124278,
+ .99999999540410731289097193313960614895889430318945,
+ .99999999885102682756267330779455410840053741619428
+ };
+
+static void fht(float *fz, int n)
+{+ int i,k1,k2,k3,k4;
+ float *fi, *fn, *gi;
+ float *tri;
+
+ fn = fi = fz + n;
+ do {+ float f0,f1,f2,f3;
+ fi -= 4;
+ f1 = fi[0]-fi[1];
+ f0 = fi[0]+fi[1];
+ f3 = fi[2]-fi[3];
+ f2 = fi[2]+fi[3];
+ fi[2] = (f0-f2);
+ fi[0] = (f0+f2);
+ fi[3] = (f1-f3);
+ fi[1] = (f1+f3);
+ } while (fi != fz);
+
+ tri = &costab[0];
+ k1 = 1;
+ do {+ float s1, c1;
+ int kx;
+ k1 *= 4;
+ k2 = k1 << 1;
+ kx = k1 >> 1;
+ k4 = k2 << 1;
+ k3 = k2 + k1;
+ fi = fz;
+ gi = fi + kx;
+ do {+ float f0,f1,f2,f3;
+ f1 = fi[0] - fi[k1];
+ f0 = fi[0] + fi[k1];
+ f3 = fi[k2] - fi[k3];
+ f2 = fi[k2] + fi[k3];
+ fi[k2] = f0 - f2;
+ fi[0 ] = f0 + f2;
+ fi[k3] = f1 - f3;
+ fi[k1] = f1 + f3;
+ f1 = gi[0] - gi[k1];
+ f0 = gi[0] + gi[k1];
+ f3 = SQRT2 * gi[k3];
+ f2 = SQRT2 * gi[k2];
+ gi[k2] = f0 - f2;
+ gi[0 ] = f0 + f2;
+ gi[k3] = f1 - f3;
+ gi[k1] = f1 + f3;
+ gi += k4;
+ fi += k4;
+ } while (fi<fn);
+ c1 = tri[0];
+ s1 = tri[1];
+ for (i = 1; i < kx; i++) {+ float c2,s2;
+ c2 = 1.0 - 2.0*s1*s1;
+ s2 = 2.0*s1*c1;
+ fi = fz + i;
+ gi = fz + k1 - i;
+ do {+ float a,b,g0,f0,f1,g1,f2,g2,f3,g3;
+ b = s2*fi[k1] - c2*gi[k1];
+ a = c2*fi[k1] + s2*gi[k1];
+ f1 = fi[0 ] - a;
+ f0 = fi[0 ] + a;
+ g1 = gi[0 ] - b;
+ g0 = gi[0 ] + b;
+ b = s2*fi[k3] - c2*gi[k3];
+ a = c2*fi[k3] + s2*gi[k3];
+ f3 = fi[k2] - a;
+ f2 = fi[k2] + a;
+ g3 = gi[k2] - b;
+ g2 = gi[k2] + b;
+ b = s1*f2 - c1*g3;
+ a = c1*f2 + s1*g3;
+ fi[k2] = f0 - a;
+ fi[0 ] = f0 + a;
+ gi[k3] = g1 - b;
+ gi[k1] = g1 + b;
+ b = c1*g2 - s1*f3;
+ a = s1*g2 + c1*f3;
+ gi[k2] = g0 - a;
+ gi[0 ] = g0 + a;
+ fi[k3] = f1 - b;
+ fi[k1] = f1 + b;
+ gi += k4;
+ fi += k4;
+ } while (fi<fn);
+ c2 = c1;
+ c1 = c2 * tri[0] - s1 * tri[1];
+ s1 = c2 * tri[1] + s1 * tri[0];
+ }
+ tri += 2;
+ } while (k4<n);
+}
+
+void fft(float *x_real, float *energy, int N)
+{+ float a,b;
+ int i,j;
+
+ fht(x_real,N);
+
+ energy[0] = x_real[0] * x_real[0];
+ for (i=1,j=N-1;i<N/2;i++,j--) {+ a = x_real[i];
+ b = x_real[j];
+ energy[i]=(a*a + b*b)/2;
+
+ if (energy[i] < 0.0005f) {+ energy[i] = 0.0005f;
+ }
+ }
+ energy[N/2] = x_real[N/2] * x_real[N/2];
+}
+
+void fft_side(float x_real0[],float x_real1[], float *x_real, float *energy, int N, int sign)
+/*
+x_real0 = x_real from channel 0
+x_real1 = x_real from channel 1
+sign = 1: compute mid channel energy, ax, bx
+sign =-1: compute side channel energy, ax, bx
+*/
+{+ float a,b;
+ int i,j;
+
+#define XREAL(j) ((x_real0[j]+sign*x_real1[j])/SQRT2)
+
+ energy[0] = XREAL(0) * XREAL(0);
+ x_real[0] = XREAL(0);
+
+ for (i=1,j=N-1;i<N/2;i++,j--) {+ a = x_real[i] = XREAL(i);
+ b = x_real[j] = XREAL(j);
+ energy[i]=(a*a + b*b)/2;
+
+ if (energy[i] < 0.0005) {+ energy[i] = 0.0005;
+ x_real[i] = x_real[j] = 0.0223606797749978970790696308768019662239; /* was sqrt(0.0005) */
+ }
+ }
+ energy[N/2] = XREAL(N/2) * XREAL(N/2);
+}
--- a/psych.c
+++ b/psych.c
@@ -52,9 +52,9 @@
Source file:
-$Id: psych.c,v 1.6 1999/12/22 14:35:15 menno Exp $
-$Id: psych.c,v 1.6 1999/12/22 14:35:15 menno Exp $
-$Id: psych.c,v 1.6 1999/12/22 14:35:15 menno Exp $
+$Id: psych.c,v 1.7 1999/12/23 14:34:51 menno Exp $
+$Id: psych.c,v 1.7 1999/12/23 14:34:51 menno Exp $
+$Id: psych.c,v 1.7 1999/12/23 14:34:51 menno Exp $
**********************************************************************/
@@ -65,6 +65,7 @@
#include "tf_main.h"
#include "psych.h"
+void fft(float *x_real, float *energy, int N);
SR_INFO sr_info_aac[MAX_SAMPLING_RATES+1] =
{@@ -268,8 +269,6 @@
double sample[MAX_TIME_CHANNELS+2][BLOCK_LEN_LONG*2];
/* sample value */
-FFT_TABLE_LONG fft_tbl_long; /* table for long fft */
-FFT_TABLE_SHORT fft_tbl_short; /* table for short fft */
PARTITION_TABLE_LONG *part_tbl_long;
/* partition table for long block */
PARTITION_TABLE_SHORT *part_tbl_short;
@@ -279,116 +278,19 @@
PSY_STATVARIABLE_SHORT psy_stvar_short[MAX_TIME_CHANNELS+2];
/* static variables for short block */
/* added by T. Araki (1997.10.16) end */
+DYN_PARTITION_TABLE_SHORT dyn_short;
+DYN_PARTITION_TABLE_LONG dyn_long;
void EncTf_psycho_acoustic_init( void )
{int chanNum;
- /* added by T. Araki (1997.10.16) */
- psy_fft_table_init(&fft_tbl_long, &fft_tbl_short);
- /* initializing fft table */
for (chanNum=0;chanNum<MAX_TIME_CHANNELS+2;chanNum++) {psy_calc_init(&sample[chanNum], &psy_stvar_long[chanNum], &psy_stvar_short[chanNum]);
/* initializing static variables */
}
- /* added by T. Araki (1997.10.16) end */
}
-/* added by T. Okada (1997.07.10) */
-void psy_fft_table_init(FFT_TABLE_LONG *fft_tbl_long,
- FFT_TABLE_SHORT *fft_tbl_short
- )
-{ -
- int i,j,k,n,n2,n4,n8;
- double c,s,dc,ds,t;
-
- /* generating Hann window */
- for(i = 0; i < BLOCK_LEN_LONG*2; i++)
- fft_tbl_long->hw[i] = 0.5 * (1-cos(2.0*M_PI*(i+0.5)/(BLOCK_LEN_LONG*2)));
- for(i = 0; i < BLOCK_LEN_SHORT*2; i++)
- fft_tbl_short->hw[i] = 0.5 * (1-cos(2.0*M_PI*(i+0.5)/(BLOCK_LEN_SHORT*2)));
-
- /* generating sin table (long) */
- n = BLOCK_LEN_LONG * 2;
- n2 = n/2;
- n4 = n2/2;
- n8 = n4/2;
-
- t = sin(M_PI / (double)n);
- dc = 2.0*t*t;
- ds = sqrt(dc * (2.0 - dc));
- t = 2*dc;
- c = fft_tbl_long->st[n4] = 1.0;
- s = fft_tbl_long->st[0] = 0;
-
- for(i = 1; i < n8; i++){- c -= dc; dc += t * c;
- s += ds; ds -= t * s;
- fft_tbl_long->st[i] = s; fft_tbl_long->st[n4 - i] = c;
- }
- if (n8 != 0) fft_tbl_long->st[n8] = sqrt(0.5);
- for (i = 0; i < n4; i++)
- fft_tbl_long->st[n2 - i] = fft_tbl_long->st[i];
- for (i = 0; i < n2 + n4; i++)
- fft_tbl_long->st[i + n2] = -fft_tbl_long->st[i];
-
- /* generating sin table (short) */
- n = BLOCK_LEN_SHORT * 2;
- n2 = n/2;
- n4 = n2/2;
- n8 = n4/2;
-
- t = sin(M_PI / (double)n);
- dc = 2*t*t;
- ds = sqrt(dc * (2.0 - dc));
- t = 2*dc;
- c = fft_tbl_short->st[n4] = 1.0;
- s = fft_tbl_short->st[0] = 0;
-
- for(i = 1; i < n8; i++){- c -= dc; dc += t * c;
- s += ds; ds -= t * s;
- fft_tbl_short->st[i] = s; fft_tbl_short->st[n4 - i] = c;
- }
- if (n8 != 0) fft_tbl_short->st[n8] = sqrt(0.5);
- for (i = 0; i < n4; i++)
- fft_tbl_short->st[n2 - i] = fft_tbl_short->st[i];
- for (i = 0; i < n2 + n4; i++)
- fft_tbl_short->st[i + n2] = - fft_tbl_short->st[i];
-
- /* generating bit inverse table (long) */
- n = BLOCK_LEN_LONG * 2;
- n2 = n/2; i = j = 0;
-
- for(;;){- fft_tbl_long->brt[i] = j;
- if( ++i >= n ) break;
- k = n2;
- while(k <= j){- j -= k;
- k /= 2;
- }
- j += k;
- }
-
- /* generating bit inverse table (short) */
- n = BLOCK_LEN_SHORT * 2;
- n2 = n/2; i = j = 0;
-
- for(;;){- fft_tbl_short->brt[i] = j;
- if( ++i >= n ) break;
- k = n2;
- while(k <= j){- j -= k;
- k /= 2;
- }
- j += k;
- }
-}
-/* added by T. Okada (1997.07.10) end */
-
/* added by T. Araki (1997.07.10) */
void psy_part_table_init(
double sampling_rate,
@@ -426,7 +328,8 @@
double ji = j + ((*part_tbl_long)->width[b]-1)/2.0;
double freq = (*part_tbl_long)->sampling_rate*ji/2048;
double bark = 13*atan(0.00076*freq)+3.5*atan((freq/7500)*(freq/7500));
- (*part_tbl_long)->bval[b] = bark;
+ //(*part_tbl_long)->bval[b] = bark;
+ dyn_long.bval[b] = bark;
}
}
@@ -442,7 +345,7 @@
double ji = j + ((*part_tbl_short)->width[b]-1)/2.0;
double freq = (*part_tbl_short)->sampling_rate*ji/256;
double bark = 13*atan(0.00076*freq) + 3.5*atan((freq/7500)*(freq/7500));
- (*part_tbl_short)->bval[b]=bark;
+ dyn_short.bval[b]=bark;
}
}
@@ -452,9 +355,11 @@
int b, bb;
for( b = 0; b < (*part_tbl_long)->len; b++) {- b2 = (*part_tbl_long)->bval[b];
+ //b2 = (*part_tbl_long)->bval[b];
+ b2 = dyn_long.bval[b];
for( bb = 0; bb < (*part_tbl_long)->len; bb++) {- b1 = (*part_tbl_long)->bval[bb];
+ //b1 = (*part_tbl_long)->bval[bb];
+ b1 = dyn_long.bval[bb];
//tmpx = (b2 >= b1) ? 3.0*(b2-b1) : 1.5*(b2-b1);
tmpx = (bb >= b) ? 3.0*(b2-b1) : 1.5*(b2-b1);
@@ -463,14 +368,14 @@
tmpy = 15.811389 + 7.5*(tmpx + 0.474)-17.5 *sqrt(1.0 + (tmpx+0.474)*(tmpx+0.474));
- (*part_tbl_long)->spreading[b][bb] = ( tmpy < -100.0 ? 0.0 : pow(10.0, (tmpz + tmpy)/10.0) );
+ dyn_long.spreading[b][bb] = ( tmpy < -100.0 ? 0.0 : pow(10.0, (tmpz + tmpy)/10.0) );
}
}
for( b = 0; b < (*part_tbl_short)->len; b++) {- b2 = (*part_tbl_short)->bval[b];
+ b2 = dyn_short.bval[b];
for( bb = 0; bb < (*part_tbl_short)->len; bb++) {- b1 = (*part_tbl_short)->bval[bb];
+ b1 = dyn_short.bval[bb];
//tmpx = (b2 >= b1) ? 3.0*(b2-b1) : 1.5*(b2-b1);
tmpx = (bb >= b) ? 3.0*(b2-b1) : 1.5*(b2-b1);
@@ -479,7 +384,7 @@
tmpy = 15.811389 + 7.5*(tmpx + 0.474)-17.5 *sqrt(1.0 + (tmpx+0.474)*(tmpx+0.474));
- (*part_tbl_short)->spreading[b][bb] = ( tmpy < -100.0 ? 0.0 : pow(10.0, (tmpz + tmpy)/10.0) );
+ dyn_short.spreading[b][bb] = ( tmpy < -100.0 ? 0.0 : pow(10.0, (tmpz + tmpy)/10.0) );
}
}
}
@@ -489,8 +394,8 @@
tmp = 0.0;
for( bb = 0; bb < (*part_tbl_long)->len; bb++)
//tmp += sprdngf( (*part_tbl_long),(*part_tbl_short), bb, b, 0);
- tmp += (*part_tbl_long)->spreading[bb][b];
- (*part_tbl_long)->rnorm[b] = 1.0/tmp;
+ tmp += dyn_long.spreading[bb][b];
+ dyn_long.rnorm[b] = 1.0/tmp;
}
/* added by T. Okada (1997.07.10) end */
@@ -499,17 +404,26 @@
tmp = 0.0;
for( bb = 0; bb < (*part_tbl_short)->len; bb++)
//tmp += sprdngf( (*part_tbl_long), (*part_tbl_short), bb, b, 1);
- tmp += (*part_tbl_short)->spreading[bb][b];
- (*part_tbl_short)->rnorm[b] = 1.0/tmp;
+ tmp += dyn_short.spreading[bb][b];
+ dyn_short.rnorm[b] = 1.0/tmp;
}
/* added by T. Araki (1997.10.16) end */
for(b = 0; b < (*part_tbl_long)->len; b++) {- (*part_tbl_long)->bmax[b] = pow(10, -3*(0.5+0.5*(M_PI*(min((*part_tbl_long)->bval[b], 15.5)/15.5))));
+ dyn_long.bmax[b] = pow(10, -3*(0.5+0.5*(M_PI*(min(dyn_long.bval[b], 15.5)/15.5))));
}
for(b = 0; b < (*part_tbl_short)->len; b++) {- (*part_tbl_short)->bmax[b] = pow(10, -3*(0.5+0.5*(M_PI*(min((*part_tbl_short)->bval[b], 15.5)/15.5))));
+ dyn_short.bmax[b] = pow(10, -3*(0.5+0.5*(M_PI*(min(dyn_short.bval[b], 15.5)/15.5))));
}
+
+ /* generating Hann window */
+ for(b = 0; b < BLOCK_LEN_LONG*2; b++)
+ dyn_long.hw[b] = 0.5 * (1-cos(2.0*M_PI*(b+0.5)/(BLOCK_LEN_LONG*2)));
+ for(b = 0; b < BLOCK_LEN_SHORT*2; b++)
+ dyn_short.hw[b] = 0.5 * (1-cos(2.0*M_PI*(b+0.5)/(BLOCK_LEN_SHORT*2)));
+
+ (*part_tbl_short)->dyn = &dyn_short;
+ (*part_tbl_long)->dyn = &dyn_long;
}
@@ -526,15 +440,6 @@
sample[ch][i] = 0.0;
}
- /* for(ch = 0; ch < Chans; ++ch){ */- for(i = 0; i < BLOCK_LEN_LONG*3; i++){- psy_stvar_long->fft_r[i] = 0.0;
- psy_stvar_long->fft_f[i] = 0.0;
- }
- /* }*/
-
- psy_stvar_long->p_fft = 0;
-
/* for(ch = 0; ch < Chans; ++ch){*/ for(i = 0; i < NPART_LONG*2; i++){psy_stvar_long->nb[i] = 90.0;
@@ -544,16 +449,6 @@
psy_stvar_long->p_nb = NPART_LONG;
/* added by T. Araki (1997.07.10) end */
-/* added by T. Araki (1997.10.16) */
- /* for(ch = 0; ch < Chans; ++ch){*/- for(i = 0; i < BLOCK_LEN_SHORT; i++) {- psy_stvar_short->last6_fft_r[i] = 0.0;
- psy_stvar_short->last6_fft_f[i] = 0.0;
- psy_stvar_short->last7_fft_r[i] = 0.0;
- psy_stvar_short->last7_fft_f[i] = 0.0;
- }
- /* }*/
-
/* for(ch = 0; ch < Chans; ++ch){*/ for(i = 0; i < NPART_SHORT; i++){psy_stvar_short->last7_nb[i] = 90.0;
@@ -625,10 +520,11 @@
{ch = 0;
psy_step1(p_time_signal,sample, no_of_chan);
- psy_step2(&sample[no_of_chan], &psy_stvar_long[no_of_chan], &psy_stvar_short[no_of_chan], &fft_tbl_long,
- &fft_tbl_short, ch);
- psy_step3(&psy_stvar_long[no_of_chan], &psy_stvar_short[no_of_chan], &psy_var_long, &psy_var_short, ch);
- psy_step4(&psy_stvar_long[no_of_chan], &psy_stvar_short[no_of_chan], &psy_var_long, &psy_var_short, ch);
+ psy_step2(&sample[no_of_chan], part_tbl_long, part_tbl_short, psy_stvar_long,
+ psy_stvar_short, ch);
+ if (no_of_chan < 2) {+ psy_step4(&psy_stvar_long[no_of_chan], &psy_stvar_short[no_of_chan], &psy_var_long, &psy_var_short, no_of_chan);
+ }
if (no_of_chan == 0) {for (b = 0; b < 70; b++)
@@ -721,235 +617,204 @@
}
}
-void psy_step2(double sample[][BLOCK_LEN_LONG*2],
- PSY_STATVARIABLE_LONG *psy_stvar_long,
+void psy_step2(double sample[][BLOCK_LEN_LONG*2],
+ PARTITION_TABLE_LONG *part_tbl_long,
+ PARTITION_TABLE_SHORT *part_tbl_short,
+ PSY_STATVARIABLE_LONG *psy_stvar_long,
PSY_STATVARIABLE_SHORT *psy_stvar_short,
- FFT_TABLE_LONG *fft_tbl_long,
- FFT_TABLE_SHORT *fft_tbl_short,
- int ch
+ int ch
)
{- int w,i,j,k,l,h,n,d,ik,k2,n4;
- double t,s,c,dx,dy;
- double *xl,*yl;
+ int i,l,n;
- /* FFT for long */
- xl = (double *)malloc( sizeof(double) * BLOCK_LEN_LONG * 2 );
- yl = (double *)malloc( sizeof(double) * BLOCK_LEN_LONG * 2 );
+ if (ch < 2) {+ for(i = 0; i < BLOCK_LEN_LONG*2; i++){+ psy_stvar_long[ch].xreal[i] = part_tbl_long->dyn->hw[i] * sample[ch][i];
+ }
- psy_stvar_long->p_fft += BLOCK_LEN_LONG;
+ n = BLOCK_LEN_LONG*2;
- if(psy_stvar_long->p_fft == BLOCK_LEN_LONG * 3)
- psy_stvar_long->p_fft = 0;
+ fft(psy_stvar_long[ch].xreal, psy_stvar_long[ch].fft_energy, n);
- /* window *//* static int it = 0; */
- for(i = 0; i < BLOCK_LEN_LONG*2; i++){- xl[i] = fft_tbl_long->hw[i] * sample[ch][i];
- yl[i] = 0.0;
- }
+ for(l = 0; l < MAX_SHORT_WINDOWS; l++){- n = BLOCK_LEN_LONG*2;
- n4 = n/4;
-
- for (i = 0; i < n; i++) { /* bit inverse */- j = fft_tbl_long->brt[i];
- if (i < j) {- t = xl[i]; xl[i] = xl[j]; xl[j] = t;
- t = yl[i]; yl[i] = yl[j]; yl[j] = t;
- }
- }
- for (k = 1; k < n; k = k2) { /* translation */- h = 0; k2 = k + k; d = n / k2;
- for (j = 0; j < k; j++) {- c = fft_tbl_long->st[h + n4];
- s = fft_tbl_long->st[h];
- for (i = j; i < n; i += k2) {- ik = i + k;
- dx = s * yl[ik] + c * xl[ik];
- dy = c * yl[ik] - s * xl[ik];
- xl[ik] = xl[i] - dx; xl[i] += dx;
- yl[ik] = yl[i] - dy; yl[i] += dy;
+ /* window */
+ for(i = 0; i < BLOCK_LEN_SHORT*2; i++){+ psy_stvar_short[ch].xreal[l][i] = part_tbl_short->dyn->hw[i] * sample[ch][/*OFFSET_FOR_SHORT +*/ BLOCK_LEN_SHORT * l + i];
}
- h += d;
- }
- }
- for(w = 0; w < BLOCK_LEN_LONG; w++){- psy_stvar_long->fft_r[w+psy_stvar_long->p_fft]
- = sqrt(xl[w]*xl[w] + yl[w]*yl[w]);
+ n = BLOCK_LEN_SHORT*2;
- if( xl[w] > 0.0 ){- if( yl[w] >= 0.0 )
- psy_stvar_long->fft_f[w+psy_stvar_long->p_fft] = atan( yl[w] / xl[w] );
- else
- psy_stvar_long->fft_f[w+psy_stvar_long->p_fft] = atan( yl[w] / xl[w] )+ M_PI * 2.0;
- } else if( xl[w] < 0.0 ) {- psy_stvar_long->fft_f[w+psy_stvar_long->p_fft] = atan( yl[w] / xl[w] ) + M_PI;
- } else {- if( yl[w] > 0.0 )
- psy_stvar_long->fft_f[w+psy_stvar_long->p_fft] = M_PI * 0.5;
- else if( yl[w] < 0.0 )
- psy_stvar_long->fft_f[w+psy_stvar_long->p_fft] = M_PI * 1.5;
- else
- psy_stvar_long->fft_f[w+psy_stvar_long->p_fft] = 0.0; /* tmp */
+ fft(psy_stvar_short[ch].xreal[l], psy_stvar_short[ch].fft_energy[l], n);
}
- }
+ } else {- if (xl) free(xl);
- if (yl) free(yl);
+ n = BLOCK_LEN_LONG*2;
+ fft_side(psy_stvar_long[0].xreal, psy_stvar_long[1].xreal, psy_stvar_long[ch].xreal,
+ psy_stvar_long[ch].fft_energy, n, (ch==2)?(1):(-1));
- /* added by T. Araki (1997.10.16) */
- /* FFT for short */
- xl = (double *)malloc( sizeof(double) * BLOCK_LEN_SHORT * 2 );
- yl = (double *)malloc( sizeof(double) * BLOCK_LEN_SHORT * 2 );
+ for(l = 0; l < MAX_SHORT_WINDOWS; l++){- for(l = 0; l < MAX_SHORT_WINDOWS; l++){-
- /* window */
- for(i = 0; i < BLOCK_LEN_SHORT*2; i++){- xl[i] = fft_tbl_short->hw[i] * sample[ch][/*OFFSET_FOR_SHORT +*/ BLOCK_LEN_SHORT * l + i];
- yl[i] = 0.0;
+ n = BLOCK_LEN_SHORT*2;
+
+ fft_side(psy_stvar_short[0].xreal[l], psy_stvar_short[1].xreal[l], psy_stvar_short[ch].xreal[l],
+ psy_stvar_short[ch].fft_energy[l], n, (ch==2)?(1):(-1));
}
+ }
+}
- n = BLOCK_LEN_SHORT*2;
- n4 = n/4;
+void psy_step4(PSY_STATVARIABLE_LONG *psy_stvar_long,
+ PSY_STATVARIABLE_SHORT *psy_stvar_short,
+ PSY_VARIABLE_LONG *psy_var_long,
+ PSY_VARIABLE_SHORT *psy_var_short,
+ int ch
+ )
+{+ int j,i;
+ static double ax_sav[4][2][8], bx_sav[4][2][8], rx_sav[4][2][8];
+ static double ax_sav_s[4][8][2][128], bx_sav_s[4][8][2][128], rx_sav_s[4][8][2][128];
+ static int init = 1;
- for (i = 0; i < n; i++) { /* bit inverse */- j = fft_tbl_short->brt[i];
- if (i < j) {- t = xl[i]; xl[i] = xl[j]; xl[j] = t;
- t = yl[i]; yl[i] = yl[j]; yl[j] = t;
+ if (init) {+ memset(rx_sav,0, sizeof(rx_sav));
+ memset(ax_sav,0, sizeof(ax_sav));
+ memset(bx_sav,0, sizeof(bx_sav));
+ memset(rx_sav_s,0, sizeof(rx_sav_s));
+ memset(ax_sav_s,0, sizeof(ax_sav_s));
+ memset(bx_sav_s,0, sizeof(bx_sav_s));
+ init = 0;
+ }
+
+ for(i = 0; i < MAX_SHORT_WINDOWS; i++){+ for(j = 0; j < BLOCK_LEN_SHORT; j++){+ double an, a1, a2;
+ double bn, b1, b2;
+ double rn, r1, r2;
+ double numre, numim, den;
+
+ a2 = ax_sav_s[ch][i][1][j];
+ b2 = bx_sav_s[ch][i][1][j];
+ r2 = rx_sav_s[ch][i][1][j];
+ a1 = ax_sav_s[ch][i][1][j] = ax_sav_s[ch][i][0][j];
+ b1 = bx_sav_s[ch][i][1][j] = bx_sav_s[ch][i][0][j];
+ r1 = rx_sav_s[ch][i][1][j] = rx_sav_s[ch][i][0][j];
+ an = ax_sav_s[ch][i][0][j] = psy_stvar_short->xreal[i][j];
+ bn = bx_sav_s[ch][i][0][j] = j==0 ? psy_stvar_short->xreal[i][0] : psy_stvar_short->xreal[i][BLOCK_LEN_SHORT-j];
+ rn = rx_sav_s[ch][i][0][j] = sqrt(psy_stvar_short->fft_energy[i][j]);
+
+ { /* square (x1,y1) */+ if( r1 != 0.0 ) {+ numre = (a1*b1);
+ numim = (a1*a1-b1*b1)*0.5;
+ den = r1*r1;
+ } else {+ numre = 1.0;
+ numim = 0.0;
+ den = 1.0;
+ }
}
- }
- for (k = 1; k < n; k = k2) { /* translation */- h = 0; k2 = k + k; d = n / k2;
- for (j = 0; j < k; j++) {- c = fft_tbl_short->st[h + n4];
- s = fft_tbl_short->st[h];
- for (i = j; i < n; i += k2) {- ik = i + k;
- dx = s * yl[ik] + c * xl[ik];
- dy = c * yl[ik] - s * xl[ik];
- xl[ik] = xl[i] - dx; xl[i] += dx;
- yl[ik] = yl[i] - dy; yl[i] += dy;
+
+ { /* multiply by (x2,-y2) */+ if( r2 != 0.0 ) {+ double tmp2 = (numim+numre)*(a2+b2)*0.5;
+ double tmp1 = -a2*numre+tmp2;
+ numre = -b2*numim+tmp2;
+ numim = tmp1;
+ den *= r2;
+ } else {+ /* do nothing */
}
- h += d;
}
- }
- for(w = 0; w < BLOCK_LEN_SHORT; w++){- psy_stvar_short->fft_r[l][w]
- = sqrt(xl[w]*xl[w] + yl[w]*yl[w]);
+ { /* r-prime factor */+ double tmp = (2.0*r1-r2)/den;
+ numre *= tmp;
+ numim *= tmp;
+ }
- if( xl[w] > 0.0 ){- if( yl[w] >= 0.0 )
- psy_stvar_short->fft_f[l][w] = atan( yl[w] / xl[w] );
- else
- psy_stvar_short->fft_f[l][w] = atan( yl[w] / xl[w] )+ M_PI * 2.0;
- } else if( xl[w] < 0.0 ) {- psy_stvar_short->fft_f[l][w] = atan( yl[w] / xl[w] ) + M_PI;
+ if( (den=rn+fabs(2.0*r1-r2)) != 0.0 ) {+ numre = (an+bn)/2.0-numre;
+ numim = (an-bn)/2.0-numim;
+ psy_var_short->c[i][j] = sqrt(numre*numre+numim*numim)/den;
} else {- if( yl[w] > 0.0 )
- psy_stvar_short->fft_f[l][w] = M_PI * 0.5;
- else if( yl[w] < 0.0 )
- psy_stvar_short->fft_f[l][w] = M_PI * 1.5;
- else
- psy_stvar_short->fft_f[l][w] = 0.0; /* tmp */
+ psy_var_short->c[i][j] = 0.0;
}
}
}
- if (xl) free(xl);
- if (yl) free(yl);
- /* added by T. Araki (1997.10.16) end */
-}
+ for(j = 0; j < 8; j++){+ double an, a1, a2;
+ double bn, b1, b2;
+ double rn, r1, r2;
+ double numre, numim, den;
+
+ a2 = ax_sav[ch][1][j];
+ b2 = bx_sav[ch][1][j];
+ r2 = rx_sav[ch][1][j];
+ a1 = ax_sav[ch][1][j] = ax_sav[ch][0][j];
+ b1 = bx_sav[ch][1][j] = bx_sav[ch][0][j];
+ r1 = rx_sav[ch][1][j] = rx_sav[ch][0][j];
+ an = ax_sav[ch][0][j] = psy_stvar_long->xreal[j];
+ bn = bx_sav[ch][0][j] = j==0 ? psy_stvar_long->xreal[0] : psy_stvar_long->xreal[BLOCK_LEN_LONG-j];
+ rn = rx_sav[ch][0][j] = sqrt(psy_stvar_long->fft_energy[j]);
+
+ { /* square (x1,y1) */+ if( r1 != 0.0 ) {+ numre = (a1*b1);
+ numim = (a1*a1-b1*b1)*0.5;
+ den = r1*r1;
+ } else {+ numre = 1.0;
+ numim = 0.0;
+ den = 1.0;
+ }
+ }
-void psy_step3(PSY_STATVARIABLE_LONG *psy_stvar_long,
- PSY_STATVARIABLE_SHORT *psy_stvar_short,
- PSY_VARIABLE_LONG *psy_var_long,
- PSY_VARIABLE_SHORT *psy_var_short,
- int ch
- )
-{- int w,i;
- int p1_l,p2_l;
+ { /* multiply by (x2,-y2) */+ if( r2 != 0.0 ) {+ double tmp2 = (numim+numre)*(a2+b2)*0.5;
+ double tmp1 = -a2*numre+tmp2;
+ numre = -b2*numim+tmp2;
+ numim = tmp1;
+ den *= r2;
+ } else {+ /* do nothing */
+ }
+ }
- p1_l = psy_stvar_long->p_fft - BLOCK_LEN_LONG;
- if( p1_l < 0 )
- p1_l = BLOCK_LEN_LONG * 2;
- p2_l = p1_l - BLOCK_LEN_LONG;
- if( p2_l < 0 )
- p2_l = BLOCK_LEN_LONG * 2;
-
- for(w = 0; w < BLOCK_LEN_LONG; w++){- psy_var_long->r_pred[w] = 2.0 * psy_stvar_long->fft_r[p1_l + w] - psy_stvar_long->fft_r[p2_l + w];
- psy_var_long->f_pred[w] = 2.0 * psy_stvar_long->fft_f[p1_l + w] - psy_stvar_long->fft_f[p2_l + w];
- }
-
- /* added by T. Araki (1997.10.16) */
- for(w = 0; w < BLOCK_LEN_SHORT; w++){- psy_var_short->r_pred[0][w] = 2.0 * psy_stvar_short->last7_fft_r[w] - psy_stvar_short->last6_fft_r[w];
- psy_var_short->f_pred[0][w] = 2.0 * psy_stvar_short->last7_fft_f[w] - psy_stvar_short->last6_fft_f[w];
- psy_var_short->r_pred[1][w] = 2.0 * psy_stvar_short->fft_r[0][w] - psy_stvar_short->last7_fft_r[w];
- psy_var_short->f_pred[1][w] = 2.0 * psy_stvar_short->fft_f[0][w] - psy_stvar_short->last7_fft_f[w];
- }
-
- for(i = 2; i < MAX_SHORT_WINDOWS; i++){- for(w = 0; w < BLOCK_LEN_SHORT; w++){- psy_var_short->r_pred[i][w] = 2.0 * psy_stvar_short->fft_r[i - 1][w] - psy_stvar_short->fft_r[i - 2][w];
- psy_var_short->f_pred[i][w] = 2.0 * psy_stvar_short->fft_f[i - 1][w] - psy_stvar_short->fft_f[i - 2][w];
+ { /* r-prime factor */+ double tmp = (2.0*r1-r2)/den;
+ numre *= tmp;
+ numim *= tmp;
}
- }
- for(w = 0; w < BLOCK_LEN_SHORT; w++){- psy_stvar_short->last6_fft_r[w] = psy_stvar_short->fft_r[6][w];
- psy_stvar_short->last6_fft_f[w] = psy_stvar_short->fft_f[6][w];
- psy_stvar_short->last7_fft_r[w] = psy_stvar_short->fft_r[7][w];
- psy_stvar_short->last7_fft_f[w] = psy_stvar_short->fft_f[7][w];
- }
- /* added by T. Araki (1997.10.16) end */
-}
-
-void psy_step4(PSY_STATVARIABLE_LONG *psy_stvar_long,
- PSY_STATVARIABLE_SHORT *psy_stvar_short,
- PSY_VARIABLE_LONG *psy_var_long,
- PSY_VARIABLE_SHORT *psy_var_short,
- int ch
- )
-{- int w,i;
- double r,f,rp,fp;
-
- for(w = 0; w < BLOCK_LEN_LONG; w++){- r = psy_stvar_long->fft_r[psy_stvar_long->p_fft+w];
- f = psy_stvar_long->fft_f[psy_stvar_long->p_fft+w];
- rp = psy_var_long->r_pred[w];
- fp = psy_var_long->f_pred[w];
-
- if( r + fabs(rp) != 0.0 )
- psy_var_long->c[w] = sqrt( psy_sqr(r*cos(f) - rp*cos(fp))
- +psy_sqr(r*sin(f) - rp*sin(fp)) )/ ( r + fabs(rp) ) ;
- else
- psy_var_long->c[w] = 0.0; /* tmp */
- }
-
- /* added by T. Araki (1997.10.16) */
- for(i = 0; i < MAX_SHORT_WINDOWS; i++){- for(w = 0; w < BLOCK_LEN_SHORT; w++){- r = psy_stvar_short->fft_r[i][w];
- f = psy_stvar_short->fft_f[i][w];
- rp = psy_var_short->r_pred[i][w];
- fp = psy_var_short->f_pred[i][w];
-
- if( r + fabs(rp) != 0.0 )
- psy_var_short->c[i][w] = sqrt( psy_sqr(r*cos(f) - rp*cos(fp))
- +psy_sqr(r*sin(f) - rp*sin(fp)) )/ ( r + fabs(rp) ) ;
- else
- psy_var_short->c[i][w] = 0.0; /* tmp */
+ if( (den=rn+fabs(2.0*r1-r2)) != 0.0 ) {+ numre = (an+bn)/2.0-numre;
+ numim = (an-bn)/2.0-numim;
+ psy_var_long->c[j] = sqrt(numre*numre+numim*numim)/den;
+ } else {+ psy_var_long->c[j] = 0.0;
}
- }
- /* added by T. Araki (1997.10.16) end */
+ }
+ for(j = 8; j < BLOCK_LEN_LONG; j+=8){+ double tmp;
+ tmp = min(psy_var_short->c[0][j/8],psy_var_short->c[1][j/8]);
+ tmp = min(psy_var_short->c[2][j/8],tmp);
+ tmp = min(psy_var_short->c[3][j/8],tmp);
+ tmp = min(psy_var_short->c[4][j/8],tmp);
+ tmp = min(psy_var_short->c[5][j/8],tmp);
+ tmp = min(psy_var_short->c[6][j/8],tmp);
+ tmp = min(psy_var_short->c[7][j/8],tmp);
+ psy_var_long->c[j] = tmp;
+ psy_var_long->c[j+1] = tmp;
+ psy_var_long->c[j+2] = tmp;
+ psy_var_long->c[j+3] = tmp;
+ psy_var_long->c[j+4] = tmp;
+ psy_var_long->c[j+5] = tmp;
+ psy_var_long->c[j+6] = tmp;
+ psy_var_long->c[j+7] = tmp;
+ }
}
void psy_step5(PARTITION_TABLE_LONG *part_tbl_long,
@@ -970,8 +835,10 @@
/* added by T. Araki (1997.10.16) */
for(w = part_tbl_long->w_low[b]; w <= part_tbl_long->w_high[b]; w++){- psy_var_long->e[b] += psy_sqr(psy_stvar_long->fft_r[psy_stvar_long->p_fft+w]);
- tmp_cb += psy_sqr(psy_stvar_long->fft_r[psy_stvar_long->p_fft+w]) * psy_var_long->c[w];
+ //psy_var_long->e[b] += psy_sqr(psy_stvar_long->fft_r[psy_stvar_long->p_fft+w]);
+ psy_var_long->e[b] += psy_stvar_long->fft_energy[w];
+ //tmp_cb += psy_sqr(psy_stvar_long->fft_r[psy_stvar_long->p_fft+w]) * psy_var_long->c[w];
+ tmp_cb += psy_stvar_long->fft_energy[w] * psy_var_long->c[w];
}
/* added by T. Araki (1997.10.16) end */
@@ -984,9 +851,11 @@
psy_var_short->e[i][b] = 0.0;
tmp_cb = 0.0;
- for(w = part_tbl_short->w_low[b]; w <= part_tbl_short->w_high[b]; w++){- psy_var_short->e[i][b] += psy_sqr(psy_stvar_short->fft_r[i][w]);
- tmp_cb += psy_sqr(psy_stvar_short->fft_r[i][w]) * psy_var_short->c[i][w];
+ for(w = part_tbl_short->w_low[b]; w <= part_tbl_short->w_high[b]; w++) {+ //psy_var_short->e[i][b] += psy_sqr(psy_stvar_short->fft_r[i][w]);
+ psy_var_short->e[i][b] += psy_stvar_short->fft_energy[i][w];
+ //tmp_cb += psy_sqr(psy_stvar_short->fft_r[i][w]) * psy_var_short->c[i][w];
+ tmp_cb += psy_stvar_short->fft_energy[i][w] * psy_var_short->c[i][w];
}
psy_var_short->c[i][b] = tmp_cb;
@@ -1012,7 +881,7 @@
ct = 0.0;
for(bb = 0; bb < part_tbl_long->len; bb++){//sprd = sprdngf(part_tbl_long, part_tbl_short, bb, b, 0);
- sprd = part_tbl_long->spreading[bb][b];
+ sprd = part_tbl_long->dyn->spreading[bb][b];
ecb += psy_var_long->e[bb] * sprd;
ct += psy_var_long->c[bb] * sprd;
}
@@ -1021,7 +890,7 @@
} else {psy_var_long->cb[b] = 0.0;
}
- psy_stvar_long->en[b] = psy_var_long->en[b] = ecb * part_tbl_long->rnorm[b];
+ psy_stvar_long->en[b] = psy_var_long->en[b] = ecb * part_tbl_long->dyn->rnorm[b];
}
/* added by T. Araki (1997.10.16) */
@@ -1031,7 +900,7 @@
ct = 0.0;
for(bb = 0; bb < part_tbl_short->len; bb++){//sprd = sprdngf(part_tbl_long, part_tbl_short, bb, b, 1);
- sprd = part_tbl_short->spreading[bb][b];
+ sprd = part_tbl_short->dyn->spreading[bb][b];
ecb += psy_var_short->e[i][bb] * sprd;
ct += psy_var_short->c[i][bb] * sprd;
}
@@ -1040,7 +909,7 @@
} else {psy_var_short->cb[i][b] = 0.0;
}
- psy_stvar_short->en[i][b] = psy_var_short->en[i][b] = ecb * part_tbl_short->rnorm[b];
+ psy_stvar_short->en[i][b] = psy_var_short->en[i][b] = ecb * part_tbl_short->dyn->rnorm[b];
}
}
/* added by T. Araki (1997.10.16) end */
@@ -1218,12 +1087,12 @@
t = 0;
if (t>1)
t = 1/t;
- tempL = max(psy_stvar_long[0].nb[p1+b]*t, min(psy_stvar_long[0].nb[p1+b], part_tbl_long->bmax[b]*psy_stvar_long[0].en[b]));
- tempR = max(psy_stvar_long[1].nb[p1+b]*t, min(psy_stvar_long[1].nb[p1+b], part_tbl_long->bmax[b]*psy_stvar_long[1].en[b]));
+ tempL = max(psy_stvar_long[0].nb[p1+b]*t, min(psy_stvar_long[0].nb[p1+b], part_tbl_long->dyn->bmax[b]*psy_stvar_long[0].en[b]));
+ tempR = max(psy_stvar_long[1].nb[p1+b]*t, min(psy_stvar_long[1].nb[p1+b], part_tbl_long->dyn->bmax[b]*psy_stvar_long[1].en[b]));
t = min(tempL,tempR);
- tempM = min(t, max(psy_stvar_long[2].nb[p1+b], min(part_tbl_long->bmax[b]*psy_stvar_long[3].en[b], psy_stvar_long[3].nb[p1+b])));
- tempS = min(t, max(psy_stvar_long[3].nb[p1+b], min(part_tbl_long->bmax[b]*psy_stvar_long[2].en[b], psy_stvar_long[2].nb[p1+b])));
+ tempM = min(t, max(psy_stvar_long[2].nb[p1+b], min(part_tbl_long->dyn->bmax[b]*psy_stvar_long[3].en[b], psy_stvar_long[3].nb[p1+b])));
+ tempS = min(t, max(psy_stvar_long[3].nb[p1+b], min(part_tbl_long->dyn->bmax[b]*psy_stvar_long[2].en[b], psy_stvar_long[2].nb[p1+b])));
if ((psy_stvar_long[0].nb[p1+b] >= 1.58*psy_stvar_long[1].nb[p1+b])&&(psy_stvar_long[1].nb[p1+b] >= 1.58*psy_stvar_long[0].nb[p1+b])) {psy_stvar_long[2].nb[p1+b] = tempM;
@@ -1241,12 +1110,12 @@
t = 0;
if (t>1)
t = 1/t;
- tempL = max(psy_stvar_short[0].nb[i][b]*t, min(psy_stvar_short[0].nb[i][b], part_tbl_short->bmax[b]*psy_stvar_short[0].en[i][b]));
- tempR = max(psy_stvar_short[1].nb[i][b]*t, min(psy_stvar_short[1].nb[i][b], part_tbl_short->bmax[b]*psy_stvar_short[1].en[i][b]));
+ tempL = max(psy_stvar_short[0].nb[i][b]*t, min(psy_stvar_short[0].nb[i][b], part_tbl_short->dyn->bmax[b]*psy_stvar_short[0].en[i][b]));
+ tempR = max(psy_stvar_short[1].nb[i][b]*t, min(psy_stvar_short[1].nb[i][b], part_tbl_short->dyn->bmax[b]*psy_stvar_short[1].en[i][b]));
t = min(tempL,tempR);
- tempM = min(t, max(psy_stvar_short[2].nb[i][b], min(part_tbl_short->bmax[b]*psy_stvar_short[3].en[i][b], psy_stvar_short[3].nb[i][b])));
- tempS = min(t, max(psy_stvar_short[3].nb[i][b], min(part_tbl_short->bmax[b]*psy_stvar_short[2].en[i][b], psy_stvar_short[2].nb[i][b])));
+ tempM = min(t, max(psy_stvar_short[2].nb[i][b], min(part_tbl_short->dyn->bmax[b]*psy_stvar_short[3].en[i][b], psy_stvar_short[3].nb[i][b])));
+ tempS = min(t, max(psy_stvar_short[3].nb[i][b], min(part_tbl_short->dyn->bmax[b]*psy_stvar_short[2].en[i][b], psy_stvar_short[2].nb[i][b])));
if ((psy_stvar_short[0].nb[i][b] >= 1.58*psy_stvar_short[1].nb[i][b])&&(psy_stvar_short[1].nb[i][b] >= 1.58*psy_stvar_short[0].nb[i][b])) {psy_stvar_short[2].nb[i][b] = tempM;
@@ -1317,7 +1186,8 @@
psy_var_long->epart[n] = 0.0;
for(w = w_low; w < w_high; w++){- psy_var_long->epart[n] += psy_sqr(psy_stvar_long->fft_r[psy_stvar_long->p_fft + w]);
+ //psy_var_long->epart[n] += psy_sqr(psy_stvar_long->fft_r[psy_stvar_long->p_fft + w]);
+ psy_var_long->epart[n] += psy_stvar_long->fft_energy[w];
}
}
@@ -1361,7 +1231,8 @@
psy_var_short->epart[i][n] = 0.0;
for(w = w_low; w < w_high; w++){- psy_var_short->epart[i][n] += psy_sqr(psy_stvar_short->fft_r[i][w]);
+ //psy_var_short->epart[i][n] += psy_sqr(psy_stvar_short->fft_r[i][w]);
+ psy_var_short->epart[i][n] += psy_stvar_short->fft_energy[i][w];
}
}
--- a/psych.h
+++ b/psych.h
@@ -63,16 +63,14 @@
/* added by T. Araki (1997.07.10) */
typedef struct {- double hw[BLOCK_LEN_LONG*2]; /* Hann window table */
- double st[BLOCK_LEN_LONG*5/2]; /* sin table*/
- int brt[BLOCK_LEN_LONG*2]; /* bit inverse table */
-} FFT_TABLE_LONG;
-
-typedef struct {- double hw[BLOCK_LEN_SHORT*2]; /* Hann window table */
- double st[BLOCK_LEN_SHORT*5/2]; /* sin table*/
- int brt[BLOCK_LEN_SHORT*2]; /* bit inverse table */
-} FFT_TABLE_SHORT;
+ double bval[NPART_LONG];
+ double qsthr[NPART_LONG];
+ double e_qsthr[NPART_LONG]; /* absolute threshold (energy) in each partition */
+ double rnorm[NPART_LONG];
+ double bmax[NPART_LONG];
+ double spreading[NPART_LONG][NPART_LONG];
+ double hw[BLOCK_LEN_LONG*2];
+} DYN_PARTITION_TABLE_LONG;
typedef struct {int sampling_rate;
@@ -80,20 +78,10 @@
int w_low[NPART_LONG];
int w_high[NPART_LONG];
int width[NPART_LONG];
- double bval[NPART_LONG];
- double qsthr[NPART_LONG];
- double e_qsthr[NPART_LONG]; /* absolute threshold (energy) in each partition */
- double rnorm[NPART_LONG];
- double bmax[NPART_LONG];
- double spreading[NPART_LONG][NPART_LONG];
+ DYN_PARTITION_TABLE_LONG *dyn;
} PARTITION_TABLE_LONG;
typedef struct {- int sampling_rate;
- int len; /* length of the table */
- int w_low[NPART_SHORT];
- int w_high[NPART_SHORT];
- int width[NPART_SHORT];
double bval[NPART_SHORT];
double qsthr[NPART_SHORT];
double e_qsthr[NPART_SHORT]; /* absolute threshold (energy) in each partition */
@@ -100,12 +88,21 @@
double rnorm[NPART_SHORT];
double bmax[NPART_SHORT];
double spreading[NPART_SHORT][NPART_SHORT];
+ double hw[BLOCK_LEN_SHORT*2];
+} DYN_PARTITION_TABLE_SHORT;
+
+typedef struct {+ int sampling_rate;
+ int len; /* length of the table */
+ int w_low[NPART_SHORT];
+ int w_high[NPART_SHORT];
+ int width[NPART_SHORT];
+ DYN_PARTITION_TABLE_SHORT *dyn;
} PARTITION_TABLE_SHORT;
typedef struct {- double fft_r[BLOCK_LEN_LONG*3];
- double fft_f[BLOCK_LEN_LONG*3];
- int p_fft; /* pointer for fft_r and fft_f */
+ float xreal[BLOCK_LEN_LONG*2];
+ float fft_energy[BLOCK_LEN_LONG*2];
double nb[NPART_LONG*2];
double en[NPART_LONG];
double save_npart_l[NSFB_LONG];
@@ -114,8 +111,6 @@
} PSY_STATVARIABLE_LONG;
typedef struct {- double r_pred[BLOCK_LEN_LONG];
- double f_pred[BLOCK_LEN_LONG];
double c[BLOCK_LEN_LONG];
double e[NPART_LONG];
double cw[NPART_LONG];
@@ -131,12 +126,8 @@
} PSY_VARIABLE_LONG;
typedef struct {- double fft_r[MAX_SHORT_WINDOWS][BLOCK_LEN_SHORT];
- double fft_f[MAX_SHORT_WINDOWS][BLOCK_LEN_SHORT];
- double last6_fft_r[BLOCK_LEN_SHORT];
- double last6_fft_f[BLOCK_LEN_SHORT];
- double last7_fft_r[BLOCK_LEN_SHORT];
- double last7_fft_f[BLOCK_LEN_SHORT];
+ float xreal[MAX_SHORT_WINDOWS][BLOCK_LEN_SHORT*2];
+ float fft_energy[MAX_SHORT_WINDOWS][BLOCK_LEN_SHORT*2];
double nb[MAX_SHORT_WINDOWS][NPART_SHORT];
double en[MAX_SHORT_WINDOWS][NPART_SHORT];
double save_npart_s[MAX_SHORT_WINDOWS][NSFB_SHORT];
@@ -145,8 +136,6 @@
} PSY_STATVARIABLE_SHORT;
typedef struct {- double r_pred[MAX_SHORT_WINDOWS][BLOCK_LEN_SHORT];
- double f_pred[MAX_SHORT_WINDOWS][BLOCK_LEN_SHORT];
double c[MAX_SHORT_WINDOWS][BLOCK_LEN_SHORT];
double e[MAX_SHORT_WINDOWS][NPART_SHORT];
double cw[MAX_SHORT_WINDOWS][NPART_SHORT];
@@ -203,10 +192,6 @@
double psy_get_absthr(double f); /* Jul 8 */
-void psy_fft_table_init(FFT_TABLE_LONG *fft_tbl_long,
- FFT_TABLE_SHORT *fft_tbl_short
- );
-
void psy_part_table_init(double sampling_rate,
PARTITION_TABLE_LONG **part_tbl_long,
PARTITION_TABLE_SHORT **part_tbl_short
@@ -223,12 +208,12 @@
);
void psy_step2(double sample[][BLOCK_LEN_LONG*2],
- PSY_STATVARIABLE_LONG *psy_stvar_long,
+ PARTITION_TABLE_LONG *part_tbl_long,
+ PARTITION_TABLE_SHORT *part_tbl_short,
+ PSY_STATVARIABLE_LONG *psy_stvar_long,
PSY_STATVARIABLE_SHORT *psy_stvar_short,
- FFT_TABLE_LONG *fft_tbl_long,
- FFT_TABLE_SHORT *fft_tbl_short,
- int ch
- );
+ int ch
+ );
void psy_step3(PSY_STATVARIABLE_LONG *psy_stvar_long,
PSY_STATVARIABLE_SHORT *psy_stvar_short,