/* Spectroscopic Toolkit version 1.92 by Pieter Suurmond, september 7, 2006. Copyright (c) 2000-2006 - Pieter Suurmond Permission is hereby granted, free of charge, to any person obtaining a copy of this software and associated documentation files (the "Software"), to deal in the Software without restriction, including without limitation the rights to use, copy, modify, merge, publish, distribute, sublicense, and/or sell copies of the Software, and to permit persons to whom the Software is furnished to do so, subject to the following conditions: The above copyright notice and this permission notice shall be included in all copies or substantial portions of the Software. Any person wishing to distribute modifications to the Software is requested to send the modifications to the original developer so that they can be incorporated into the canonical version. THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE. */ #include /* Standard headers. */ #include #include #include "ST_wavelets.h" /* For addWAVELET(). */ #include "ST_pconv.h" #include "c00.h" typedef struct { double Emin, /* Energy-level statistics per element. */ Emax, Ediff; } c04_ElementStat; #define kc04_firstE (1) /* First element (for running tests, to skip elements). */ #define kc04_lastE (1) /* Last element to include in the run. */ /* For static array-allocation & stats (don't change for running-tests).*/ #define kc04_elems (96) /*--------------------------------------------------------------------------------------------*/ short c04_Regress(ENGINEp E, FILE* msg) /* ENERGY/TIME-SMEAR, START WITH HYDROGEN??? */ { c04_ElementStat ESTAT[kc04_elems + 1]; /* +1 because we index from 1 instead of 0!!! */ double EAMPS[kc04_elems + 1] = { /* Linear amplitudes (a little "sound design"): */ /*------------------------------------------------------------------------*/ /* 1..8: (not:) H He Li Be B C N O */ 0.00, 1.24, 1.00, 0.92, 0.94, 1.00, 0.60, 0.68, 0.70, /* 9..16: F Ne Na Mg Al Si P S */ 0.96, 0.70, 0.95, 0.90, 0.97, 0.60, 0.90, 0.90, /* Above amplitudes (1 to 16 inclusive) are ok, no over/underflow. */ /* With A-factor=1.0e-15 somewhere below, and with kc04_elems = 96, */ /* and we must use zero-phase aligned cosines (not sines) at TOP. */ /*------------------------------------------------------------------------*/ /* 17..24: Cl Ar K Ca Sc Ti V Cr */ 0.72, 0.56, 0.78, 0.80, 0.80, 0.78, 0.78, 0.80, /* 0.86, 0.54 OVERFLOWS @ 67.44 secs. */ /*-------------------------------------------------------------------*/ /* 25..32: Mn Fe Co Ni Cu Zn Ga Ge */ 0.50, 0.48, 0.50, 0.50, 0.54, 0.80, 0.80, 0.80, /* to be tested on overflows / peak-values... */ /* 25-32 is a heavy job! Logfile is 117 MB, be very very patient. */ /* 0.50, 0.48 Gives Overflow @ Fe @ 1.16 minutes = 69.6 secs. */ /*-------------------------------------------------------------------*/ /* 33..40: As Se Br Kr Rb Sr Y Zr */ 1.00, 1.00, 1.00, 1.00, 1.00, 1.00, 1.00, 1.00, /*-------------------------------------------------------------------*/ /* 41..48: Nb Mo Tc Ru Rh Pd Ag Cd */ 1.00, 1.00, 1.00, 1.00, 1.00, 1.00, 1.00, 1.00, /*-------------------------------------------------------------------*/ /* 49..56: In Sn Sb Te I Xe Cs Ba */ 1.00, 1.00, 1.00, 1.00, 1.00, 1.00, 1.00, 1.00, /* 57..64: La Ce Pr Nd Pm Sm Eu Gd */ 1.00, 1.00, 1.00, 1.00, 1.00, 1.00, 1.00, 1.00, /* 65..72: Tb Dy Ho Er Tm Yb Lu Hf */ 1.00, 1.00, 1.00, 1.00, 1.00, 1.00, 1.00, 1.00, /* 73..80: Ta W Re Os Ir Pt Au Hg */ 1.00, 1.00, 1.00, 1.00, 1.00, 1.00, 1.00, 1.00, /* 81..88: Tl Pb Bi Po At Rn Fr Ra */ 1.00, 1.00, 1.00, 1.00, 1.00, 1.00, 1.00, 1.00, /* 89..96: Ac Th Pa U Np Pu Am Cm */ 1.00, 1.00, 1.00, 1.00, 1.00, 1.00, 1.00, 1.00 }; FILE *pbinDatabaseFP; pieterPacked pbinBuff; derivedData more; unsigned long linesRead = 0, linesUsed = 0, linesSkipped = 0; short e = 0; int atomicNumber, ionizationLevel, swapped = 0; double audioHertz, a, dur, topStart, top, d, durA, durB; double L, Lp, Lq; /* L,Lp,Lq for logarithmic scaling. */ /*-------------------------------------------------------- Frequency & time transpositions: */ double dayTime = 24.0 * 3600.0 / 4096.0; /* Secs per day, transposed 12 octs. */ double transpositionRatio = pow(2.0, -40.0); /* 40 octaves down. */ /*-------------------------------------------------------- Smoother frequency windowing: */ double lowestFreq = 20.0; /* (trapeziod instead of brickwall) */ double aboveLowestFreq = 1.25 * lowestFreq; /* Rise to 1 within 1 major third. */ double aboveHertz = aboveLowestFreq - lowestFreq; double nyquist = 0.5*(double)getSamplerateENGINE(E); double belowNyquist = nyquist / 1.2; /* Fall to 0 within 1 minor third. */ double belowHertz = nyquist - belowNyquist; /*-------------------------------------------------------- GLOBAL freqency/amplitude stats: */ double Fmin = 1.0e31, Fmax = 0.0, /* Audio-Hertz! */ Amin = 1.0e31, Amax = 0.0; /*-------------------------------------------------------- LOCAL energy stats (per element): */ for (atomicNumber=1; atomicNumber<=kc04_elems; atomicNumber++) { /* Init statistics-array for ALL. */ ESTAT[atomicNumber].Emax = 0.0; /* Ediffs will be calculated at the */ ESTAT[atomicNumber].Emin = 1.0e31; /* first (analysis-)run. */ } /*---------------------------------------------------- Open binary version of Kurucz's file. */ pbinDatabaseFP = fopen(kST_bin_base, "rb"); if (pbinDatabaseFP) { /*------------------------------ CYCLE #1 (analysis) ------------------------------------*/ printf("'c04_Regress()' analyzing database '%s'...\n", kST_bin_base); /* FOR ALL ELEMS*/ if (msg) { fprintf(msg, "kc04_firstE = %d\n", kc04_firstE); fprintf(msg, " kc04_lastE = %d\n", kc04_lastE); fprintf(msg, "Linear amplitude-weights for these elements (ions included):\n"); for (atomicNumber=kc04_firstE; atomicNumber<=kc04_lastE; atomicNumber++) { fprintf(msg, "%3d=%.2f ", atomicNumber, EAMPS[atomicNumber]); if (!(atomicNumber & 7)) fprintf(msg, "\n"); } fprintf(msg, "\n"); } while (fread(&pbinBuff, sizeof(pbinBuff), 1, pbinDatabaseFP)) /* NOT ONLY THE SE- */ { /* LECTED ONES !!! */ atomicNumber = pbinBuff.ELEM / 100; ionizationLevel = pbinBuff.ELEM - (100 * atomicNumber); if (atomicNumber <= kc04_elems) /* ALL ATOMS & IONS UP TO kc04_elems. Do NOT use */ { /* kc04_firstE & kc04_lastE during F/A-analysis. */ moreLineData(&pbinBuff, &more); /* Calc some more data. */ audioHertz = more.frequency * transpositionRatio; /* Exclude out-of-range-freqs from analysis. */ /* There could slip through some freqs that are too weak! */ if ((audioHertz > lowestFreq) && (audioHertz < nyquist)) { if (audioHertz < Fmin) /* Audio-Hertz min & max for ALL elements & ions included. */ Fmin = audioHertz; if (audioHertz > Fmax) Fmax = audioHertz; if (more.Avalue < Amin) /* A-values (not audio yet). */ Amin = more.Avalue; if (more.Avalue > Amax) Amax = more.Avalue; if ((atomicNumber >= kc04_firstE) && (atomicNumber <= kc04_lastE)) /* Skip unused ones for E-scaling. */ { if (pbinBuff.LOW.E < ESTAT[atomicNumber].Emin) /* Per element: */ ESTAT[atomicNumber].Emin = pbinBuff.LOW.E; if (pbinBuff.LOW.E > ESTAT[atomicNumber].Emax) ESTAT[atomicNumber].Emax = pbinBuff.LOW.E; if (pbinBuff.UPP.E < ESTAT[atomicNumber].Emin) ESTAT[atomicNumber].Emin = pbinBuff.UPP.E; if (pbinBuff.UPP.E > ESTAT[atomicNumber].Emax) ESTAT[atomicNumber].Emax = pbinBuff.UPP.E; if ((pbinBuff.LOW.E >= pbinBuff.UPP.E) && (pbinBuff.ELEM != swapped)) { swapped = pbinBuff.ELEM; /* Print at least SOME swapped levels. */ if (msg) fprintf(msg, "Energy levels swapped (element %d)!\n", pbinBuff.ELEM); } } } } } /*--------------------- Print statistics: ------------*/ if (msg) fprintf(msg, "\ F and A analysis include ALL elements and ions up to %d (incl):\n\ Fmin = %.2f; Fmax = %.2f.\n\ Amin = %.8f; Amax = %.2f.\n\ E analysis PER element:\n", kc04_elems, Fmin, Fmax, Amin, Amax); for (atomicNumber=kc04_firstE; atomicNumber<=kc04_lastE; atomicNumber++) /* But print ONLY the selected ones. */ { ESTAT[atomicNumber].Ediff = ESTAT[atomicNumber].Emax - ESTAT[atomicNumber].Emin; /* Calc now for faster addWAVELET-calling. */ if (ESTAT[atomicNumber].Ediff < 0) { if (msg) fprintf(msg, "Energy-stats scaling problem element=%d!\n", atomicNumber); } if (msg) fprintf(msg, "\ element=%d: Emin=%.6f Emax=%.2f Ediff=%.2f\n", atomicNumber, ESTAT[atomicNumber].Emin, ESTAT[atomicNumber].Emax, ESTAT[atomicNumber].Ediff); } /*---------------- Linear duration compensation function: ------------------------*/ /* d(loF) = 1 1 = a loF + b a = 0.5 / (loF - hiF) */ /* d(hiF) = 0.5 0.5 = a hiF + b - 0.5 = (0.5 / (loF - hiF)) hiF + b */ /* d(F) = a F + b ------------------- b = 0.5 - (0.5 hiF / (loF - hiF)) */ /* 0.5 = a (loF - hiF) */ d = Fmin - Fmax; durA = 0.5 / d; /* A=slope and B=offset. */ durB = 0.5 - (0.5 * Fmax / d); /*---------------- Linear bass-boost-filter-function (per element): -----------------*/ /* SAME "d" as with duration may be used: highest freq (per element) gets -6.021 dB. */ /*---------------- Logarithmic duration compensation function: -----------------*/ /* 1 = Lp + Lq * log(loF) */ /* 0.5 = Lp + Lq * log(hiF) - L(F) = Lp + Lq * log(F) */ /* ------------------------------ */ /* 0.5 = Lq * log(loF / hiF) 1 = Lp + Lq * log(loF) */ /* 1 = Lp + 0.5 * log(loF) / log(loF / hiF) */ /* Lq = 0.5 / log(loF / hiF) Lp = 1 - 0.5 * log(loF) / log(loF / hiF) */ /* */ L = log(Fmin / Fmax); /* Example: When loF = 100Hz and hiF = 10kHz. */ Lq = 0.5 / L; /* Then Lp = 1.5 and Lq = -0.25. */ Lp = 1.0 - (0.5 * log(Fmin) / L); /* So L(F) = 1.5 - 0.25 * log(F). */ /*--------------------------------- Thus L(100)=1, L(10000)=0.5, L(1000)=0.75. */ if (msg) { fprintf(msg, "Using lin & log frequency-scalings:\n"); fprintf(msg, " linear: d(F) = %.6f + %6f * F\n", durB, durA); fprintf(msg, " logarithmic: L(F) = %.6f + %6f * log(F)\n", Lp, Lq); } /*-------------------- CYCLE #2: Post (audio-)wavelets into linked list. -------*/ rewind(pbinDatabaseFP); printf("'c04_Regress()' posting wavelets.\n"); if (msg) { fprintf(msg, "1/transpositionRatio = %.1f.\n", 1.0 / transpositionRatio); fprintf(msg, "dayTime = %.6f s.\n", dayTime); } while (fread(&pbinBuff, sizeof(pbinBuff), 1, pbinDatabaseFP)) { /*---------------------------------------------------- Get Atomic number and ionization level. */ atomicNumber = pbinBuff.ELEM / 100; ionizationLevel = pbinBuff.ELEM - (100 * atomicNumber); /*---------------------------------------------------- ELEMENT / ION - selector: */ if ((atomicNumber <= kc04_lastE) && (atomicNumber >= kc04_firstE) && (!ionizationLevel)) /* NO IONS FOR NOW... */ { moreLineData(&pbinBuff, &more); /* Calc some more data. */ audioHertz = more.frequency * transpositionRatio; if ((audioHertz > lowestFreq) && (audioHertz < nyquist)) { /*------------------ Less brickwall-shaped filtering: --------------*/ if (audioHertz < aboveLowestFreq) a = (audioHertz - lowestFreq) / aboveHertz; /* Rise to 1 within 1 minor third. */ else if (audioHertz > belowNyquist) a = (nyquist - audioHertz) / belowHertz; /* Fall to 0 within 1 major third. */ else a = 1.0; /*---------------- Frequency-dependant Compensations ---------------------------*/ /* Only lowest freq of one metal (section) gets unit-duration and unit-amplit., */ /* all higher frequencies sound shorter. loFreq -> 1, hiFreq -> 0.5. */ d = (durA * audioHertz) + durB; L = Lp + Lq * log(audioHertz); /* And the logarithmic variant: */ /*------------------ Bass boost-filter: -------------------------*/ a *= d * L * L; /* Lowest freq 0 dB; Highest frequency -18 dB. */ /* linear as well as logarithmic!! */ /* ----------------- Amplitude: ------------------------------------*/ /* Compress amplitude-range (per section) by approx 5th-power-root. */ a *= EAMPS[atomicNumber] * pow(1.0e-15 * more.Avalue, 0.2); /* (0.2= quite compressed) */ a /= (double)(1 + ionizationLevel); /* Ions weaker (-6dB for +1) */ /*-------------------- Amplitude treshold (23 bit) ------------------*/ if (a >= 1.192092896E-7) /* (= 1 / 2^23) */ { /*-------------------- Duration: ---------------------------------------------*/ dur = dayTime * d * pow(L, 5.0); /* Lowest freq (of all!) gets unit-duration, highest 1/64. */ dur *= pow(3.00e-7 * more.Avalue, 0.2); /* Extremely compressed, minimal influence...? */ dur /= (double)(1 + ionizationLevel); /* Ions shorter (half for +1) */ /*---------------------------------------------------------------------------*/ topStart = 27.0 + dayTime * (double)(atomicNumber - (kc04_firstE - 1)); /* (Skip silence) */ /* Lower and upper J numbers determine the starttimes of exponential decaying envelopes. In case J numbers are identical, only one wavelet is created in the middle. Otherwise, 2 wavelets are created: one at left channel, another one at the right channel. */ top = pbinBuff.LOW.J; if (pbinBuff.LOW.J == pbinBuff.UPP.J) /* Single wavelet center-panned. */ { if (addWAVELET(E, msg, 0, audioHertz, a, 0.0, top, dur, kENV_ExpDecay, 0)) /* FRQ, AMP, PAN, TOP, DUR, ENV, 0 = zero-phase sine at TOP. */ { e = 300; goto clFile; } } else /* Two wavelets: one left-panned, one right-panned. */ { if (addWAVELET(E, msg, 0, audioHertz, a, -0.6, top, dur, kENV_ExpDecay, 0)) /* FRQ, AMP, PAN, TOP, DUR, ENV, 0 = zero-phase sine at TOP. */ { e = 301; goto clFile; } top = pbinBuff.UPP.J; if (addWAVELET(E, msg, 0, audioHertz, a, 0.6, top, dur, kENV_ExpDecay, 0)) /* FRQ, AMP, PAN, TOP, DUR, ENV, 0 = zero-phase sine at TOP. */ { e = 302; goto clFile; } } linesUsed++; } else { linesSkipped++; if (msg) fprintf(msg, "SKIPPED-AMP %d: Hz=%12.6f lggf=%5.2f.\n", atomicNumber, audioHertz, pbinBuff.LOGGF); } } else { linesSkipped++; if (msg) fprintf(msg, "SKIPPED-FRQ %d: Hz=%12.6f lggf=%5.2f.\n", atomicNumber, audioHertz, pbinBuff.LOGGF); } printLineEssential(&pbinBuff, msg); /* stdout or logFP. */ printLineDerived(&more, msg); } linesRead++; } if (msg) fprintf(msg, "lines used = %ld; lines skipped = %ld.\n", linesUsed, linesSkipped); if (linesRead != kKuruczLinesTotal) { if (msg) fprintf(msg, "Error reading binary database lines!\n"); e = 2; } clFile: fclose(pbinDatabaseFP); } else { if (msg) fprintf(msg, "Cannot open binary database-file '%s'!\n", kST_bin_base); e = 1; } return e; }