' ********************************************************* ' * gmt library containing independent functions * ' * g_indep.dll * ' * public declarations are in the file g_indep.bi * ' * GMT Version 11.09 * ' * dr.Godfried-Willem Raes * ' ********************************************************* ' 31.08.2008: changed to compile with PB 9.0 : point variable name changed to ppoint ' 13.06.2009: more spectral transforms added. ' 14.06.2009: Start adding function for frequency determination in data arrays (zero cross counters) ' 06.07.2009: sort functions for spectrum analysis added. ' 07.07.2009: dft procedures modified to return half-size spektrum arrays. ' 14.07.2009: array sorting functions added. ' 30.01.2010-04.02.2010: functions for frustum calculation added ' 13.05.2010: Functions for spectrum to octave-bands added. ' 02.06.2010: C convention applied for arithmetic (i = i + x becomes i += x) ' 04.06.2010: Cepstrum coding added for namuda ' 30.03.2011: adapted to PBWin10 ' 14.10.2011: RangeLim function added. ' 30.03.2012: Functions for 1/f distributions added. ' 31.03.2012: ggd function added (largest common denominator) ' 08.04.2012: histogram functions for 1/f distributions added. ' 01.12.2012: DECR0 function added ' 17.09.2014: Formule voor veelhoeksgetallenreeksen toegevoegd ' 30.12.2015: Chaos theory formula added. ' 30.07.2016: Funktie voor differentiatie en integratie toe te voegen: ' dArray : vervangt elke waarde in een array door het verschil met zijn opvolger ' fArray : vervangt elke waarde in een array door het gemiddelde met zijn opvolger ' In the works... ' 14.03.2017: japanese format date-time conversion function added: Datim ' 30.05.2017: Generic PID function added ' 29.09.2018: Function RatNumPat added, to extract repeating patterns in rational numbers ' 30.09.2018: Function Phi() added to calculation the Euler totient of a number ' 08.02.2022: RC_cal funktie toegevoegd voor de ontlaadtijd van een RC schakeling. #COMPILER PBWIN 10 #COMPILE DLL "g_indep.dll" #OPTION VERSION5 ' compile for Windows Vista, 2000 and/or NT5 #REGISTER ALL ' This metastatement should appear only once. #DIM ALL #INCLUDE "..\winapi\g_win.inc" 'subset of Win32api.inc" #INCLUDE "g_kons.bi" ' only integer constants and strings can be declared in PB #INCLUDE "g_type.bi" ' This declares all our own structures, user defined types '#INCLUDE "g_indep.bi" ' declarations for all funcs and procs that are exported. GLOBAL WarnTime AS DWORD ' global for duration of warning windows. GLOBAL WarnText AS ASCIIZ * 8000 '511 ' test. GLOBAL La AS SINGLE ' value set on init of dll GLOBAL GrondDo AS SINGLE ' calculated on initialisation GLOBAL Domq AS SINGLE ' idem. GLOBAL LogFlag AS LONG DECLARE CALLBACK FUNCTION WarnProc AS LONG ' not exported ' declares for functions called internaly as well as external: ' the .bi file with the references to the dll lib. do not suffice any longer as declares in PBWIN8 DECLARE FUNCTION CreateDefaultIniFile (f AS STRING) AS LONG DECLARE SUB Warning (m AS STRING, OPT BYVAL t AS DWORD) DECLARE SUB P2C (P AS Polar, C AS Complex) DECLARE SUB C2P (C AS Complex, P AS Polar) DECLARE FUNCTION MagCom# (G1 AS Complex) DECLARE FUNCTION AngCom# (G1 AS Complex) DECLARE SUB MulPol (G1 AS Polar, G2 AS Polar, Gr AS Polar) DECLARE SUB DivPol (G1 AS Polar, G2 AS Polar, Gr AS Polar) DECLARE FUNCTION GetWindir () AS STRING ' 06.06.2002 DECLARE FUNCTION ArcCos (BYVAL SINGLE) AS SINGLE ' inverse of cosine DECLARE FUNCTION ArcSin (BYVAL SINGLE) AS SINGLE ' inverse of sine DECLARE FUNCTION LTrimZero (a AS STRING) AS STRING DECLARE FUNCTION Bentham (BYVAL n AS SINGLE) AS SINGLE ' 12.05.2005 DECLARE FUNCTION FrustumFlow (BYVAL traject AS SINGLE, BYVAL diam AS SINGLE, BYVAL angle AS SINGLE) AS SINGLE DECLARE FUNCTION LogFile(text$, OPT BYVAL filenam$) AS LONG ' intelligent functions for : DECLARE FUNCTION Sire_Velo2MidiNoot_old (BYVAL velo AS INTEGER) AS SINGLE DECLARE FUNCTION Sire_MidiNoot2Velo_old (BYVAL noot AS SINGLE) AS WORD DECLARE FUNCTION GetNoteNrFromKloktype(BYVAL pVaccaNotes AS kloktype PTR, lb AS DWORD, ub AS DWORD, BYVAL note AS SINGLE, BYVAL tol AS SINGLE) AS BYTE DECLARE FUNCTION WaveFreq (Samp() AS SINGLE, BYVAL samplingrate AS LONG, BYVAL noise AS SINGLE) AS SINGLE DECLARE FUNCTION WaveFreq_Dbl (Samp() AS DOUBLE, BYVAL samplingrate AS LONG, BYVAL noise AS DOUBLE) AS SINGLE DECLARE FUNCTION WaveFreq_Int (Samp() AS INTEGER, BYVAL samplingrate AS LONG, BYVAL noise AS WORD) AS SINGLE DECLARE FUNCTION WaveFreq_Lng (Samp() AS LONG, BYVAL samplingrate AS LONG, BYVAL noise AS DWORD) AS SINGLE 'DECLARE FUNCTION Sire_Play (BYVAL noot AS INTEGER, BYVAL velo AS BYTE) AS LONG FUNCTION DLLMAIN(BYVAL hInstance AS LONG, _ BYVAL fwdReason AS LONG, _ BYVAL lpvReserved AS LONG) AS LONG DLLMAIN = 1 'success DLLMain = 0 is failure SELECT CASE fwdReason CASE %DLL_PROCESS_ATTACH La = 440.00 GrondDo = La * (2!^(3!/12!))/64! ' Domq = 7.94305: ' quartertone lower for semitone bands in FFT's Domq = La * (2!^(5!/24!))/64! ' for convolutions and spectrum analysis: (these need the diapason!) CASE %DLL_PROCESS_DETACH, %DLL_THREAD_ATTACH, %DLL_THREAD_DETACH CASE ELSE DLLMAIN = %False END SELECT END FUNCTION ' ***************************************************************************** ' procedure code starts here: ' ***************************************************************************** FUNCTION CheckWinVersion () EXPORT AS LONG ' checks the windows version installed on this machine 'cfr http://msdn.microsoft.com/en-us/library/ms724834(VS.85).aspx LOCAL osv AS OSVersionInfo osv.dwOSVersionInfoSize = SIZEOF(osv) GetVersionEx osv SELECT CASE osv.dwPlatformId CASE %VER_PLATFORM_WIN32s '"Win32s on Windows 3.1" FUNCTION = 31 CASE %VER_PLATFORM_WIN32_WINDOWS SELECT CASE osv.dwminorversion CASE %False '"Windows 95 system" FUNCTION = 1995 CASE <90 '"Windows 98 system" '"Build="; LOWRD(osv.dwbuildNumber) 'FUNCTION = LOWRD(osv.dwbuildNumber) 'was >> 1999 on putty FUNCTION = 1998 'changed 010807 CASE 90 '"Windows Millennium system" '"Build="; LOWRD(osv.dwbuildNumber) = 3000 FUNCTION = 1000 END SELECT 'CASE %VER_PLATFORM_WIN32_CE ' FUNCTION = 500 ' Windows CE system CASE %VER_PLATFORM_WIN32_NT '"32 bit Windows NT" SELECT CASE osv.dwmajorversion CASE 3 ' NT3 FUNCTION = 2000 + 300 + osv.dwminorversion ' 2351 for NT 3.51 CASE 4 ' NT4 FUNCTION = 2000 + 400 + osv.dwminorversion ' 2400 for NT4.0 CASE 5 ' Windows 2000 SELECT CASE osv.dwminorversion CASE 0: FUNCTION = 2000 'win 2000 CASE 1: FUNCTION = 2001 'Xp - 2001 van het release jaar.. CASE 2: FUNCTION = 2003 'windows server 2003 END SELECT CASE 6 SELECT CASE osv.dwminorversion CASE 0: FUNCTION = 2008 'server 2008 of vista - we moeten de OSversionInfoEx gebruiken als we het onderscheid willen maken.. CASE 1: FUNCTION = 7 'server 2008r2 of win 7 END SELECT CASE ELSE FUNCTION = osv.dwmajorversion + osv.dwminorversion ' 7 for win7 -> dit was per ongeluk: maj 6, min1 - niet te onderscheiden van maj 5 min 2 dus.. END SELECT END SELECT END FUNCTION FUNCTION IsNT () EXPORT AS BYTE ' returns true on Windows 2000, Windows XP , NT FUNCTION = %False IF CheckWinVersion > 1999 THEN FUNCTION = %True END FUNCTION FUNCTION ExistFile (f AS STRING) EXPORT AS BYTE LOCAL buffer AS OFSTRUCT LOCAL filehandle AS LONG STATIC nam AS ASCIIZ * 80 STATIC oldret AS BYTE IF f = "" THEN FUNCTION = %False : EXIT FUNCTION IF f = nam THEN IF oldret THEN FUNCTION = oldret END IF END IF nam = f FUNCTION = %False 'filehandle = OpenFile (f$, buffer, %OF_EXIST OR %OF_PARSE OR %OF_PROMPT) 'filehandle = OpenFile (f, buffer, %OF_PROMPT OR %OF_EXIST) ' creates a message box if file was not found. filehandle = OPENFILE (nam, buffer, %OF_EXIST) ' returns 1 for successfull opening. IF filehandle > 0 THEN '> 0 : kl 010802 oldret = %True FUNCTION = %True END IF END FUNCTION FUNCTION IniFileName () EXPORT AS STRING ' 27.11.2000: This function returns the machine specific initialisation file used to ' start up GMT applications. STATIC tog AS DWORD LOCAL nr AS LONG STATIC filenaam AS STRING LOCAL cnt AS WORD ' LOCAL ctrlwrd% IF ISFALSE tog THEN LOCAL mpName AS ASCIIZ * %MAX_COMPUTERNAME_LENGTH + 1 'constant not found in PB's winapi file...(added) GetComputerName mpName,SIZEOF(mpName) ' WinApi function mpname = MCASE$(TRIM$(mpname)) tog = %True filenaam = TRIM$(mpname) & "_" & LCASE$($INI) IF ISFALSE ExistFile(filenaam) THEN IF ExistFile ($INI) THEN FILECOPY $INI,filenaam tog = %True ELSE IF ISFALSE cnt THEN ' in this case we have to create a default ini file... nr = CreateDefaultIniFile (filenaam) Warning "<" & filenaam & "> default ini file created>", 5000 ' MSGBOX "<" & filenaam & "> default ini file created>",,FUNCNAME$ ' ctrlwrd% = &B0001001100111111 ' ! fldcw ctrlwrd% ; guarantee preservation of extended precision - windows bug. cnt = %True END IF END IF END IF FUNCTION = filenaam ELSE FUNCTION = filenaam END IF END FUNCTION SUB DeleteIfEmpty (f AS STRING) EXPORT LOCAL retval AS LONG LOCAL nr AS LONG LOCAL buffer AS OFSTRUCT retval = OPENFILE ((f), buffer,%OF_EXIST) IF retval = %True THEN ' de file bestaat... retval = FREEFILE OPEN f FOR BINARY AS #retval nr = LOF(retval) IF nr <= %False THEN CLOSE #retval KILL f ELSE CLOSE #retval END IF END IF END SUB FUNCTION MakeLong (a AS INTEGER, b AS INTEGER) EXPORT AS LONG ' this duplicates a standard C++ Win32 macro from microsoft. LOCAL tmp AS LONG tmp = b SHIFT LEFT tmp, 16 FUNCTION = tmp OR a END FUNCTION FUNCTION B2S$ (BYVAL b%) EXPORT ' returns a 2-nibble hex-string for an integer, with leading 0 ! ' in PBWIN 8.0 this is the default behaviour... FUNCTION = HEX$(b%,2) ' b% = b% AND &H0FF ' IF b% < 16 THEN ' FUNCTION = "0" + HEX$(b%) ' ELSE ' FUNCTION = HEX$(b%) ' END IF END FUNCTION FUNCTION I2S$ (BYVAL b AS INTEGER) EXPORT ' returns a 4-nibble hex-string, with leading zero's (so always 4 characters long) FUNCTION = HEX$(b%,4) ' LOCAL s AS STRING ' s = TRIM$(HEX$(b)) ' SELECT CASE LEN(s) ' CASE 1 ' s = "000" & s ' CASE 2 ' s = "00" & s ' CASE 3 ' s = "0" & s ' END SELECT ' FUNCTION = s END FUNCTION FUNCTION mmioFOURCC (a$,b$,c$,d$) EXPORT AS DWORD LOCAL FourCC AS DWORD FourCC = ASC(d$) SHIFT LEFT FourCC, 8 FourCC = FourCC OR ASC(c$) SHIFT LEFT FourCC, 8 FourCC = FourCC OR ASC(b$) SHIFT LEFT FourCC, 8 FUNCTION= FourCC OR ASC(a$) END FUNCTION SUB SolveEqDeg2 (BYVAL a!, BYVAL b!, BYVAL c!, rts AS Roots) EXPORT 'solves equations of the type ax^2 + bx + c = 0 LOCAL discriminant# discriminant# = (b!^2!) - (4!*(a! * c!)) SELECT CASE discriminant# CASE > 0 rts.rt1.real = (-b! - (discriminant#^0.5#))/ (2# * a!) rts.rt2.real = (-b! + (discriminant#^0.5#))/ (2# * a!) rts.rt1.imag = 0 rts.rt2.imag = 0 CASE ELSE rts.rt1.real = -b! / (2# * a!) rts.rt2.real = -b!/ (2# * a!) rts.rt1.imag = -(ABS(discriminant#)^0.5#) / (2# * a!) rts.rt2.imag = (ABS(discriminant#)^0.5#) / (2# * a!) END SELECT END SUB FUNCTION TriangleNumber! (BYVAL i AS SINGLE) EXPORT FUNCTION = ((i * i) - i) / 2! END FUNCTION FUNCTION TriangleNumberLimit (BYVAL x AS SINGLE) EXPORT AS SINGLE ' this function returns the highest value for i that can be passed to TriangleNumber!(i) such that ' the number returned will not exceed x. ' thus, we have to solve the equation : x = ((i*i)-i) / 2, or, i^2 - i - 2x = 0 ' The discriminant of this equation is: d = (-1)^2 - 4*(1 * -2x) ' or: d = 1 + 8x = 8x + 1 ' thus we have only real and positive solutions: LOCAL d AS SINGLE ' LOCAL root1 AS SINGLE ' LOCAL root2 AS SINGLE d = (x * 8!) + 1! ' root1 = (1 -(d^0.5)) / 2! ' root2 = (1 +(d^0.5)) / 2! FUNCTION = (1 + (SQR(d))) /2! END FUNCTION FUNCTION PolygonNumber (BYVAL nrsides AS DWORD, BYVAL n AS DWORD) EXPORT AS DWORD ' 17.09.2014 added by gwr. for the octogon production ' for triangle numbers, nrsides = 3, for octogonal numbers nrsides = 8 etc... FUNCTION = (((nrsides/2)-1)* (n^2)) + (( 2 - (nrsides/2)) * n) END FUNCTION FUNCTION BinomialComb (BYVAL n AS WORD, BYVAL k AS WORD) EXPORT AS DWORD ' n = total number of elements ' k = number of elements in a subset of n LOCAL noemer AS DWORD LOCAL teller AS DWORD LOCAL i AS WORD IF k > n THEN Warning "function called with subsets larger than total: " + FUNCNAME$, 10000 ' i = &B0001001100111111 ' ! fldcw i ; guarantee preservation of extended precision - windows bug. FUNCTION = %False EXIT FUNCTION END IF FOR i = 1 TO k noemer *= i ' berekenen k! NEXT i teller = n * (n-1) * (n-2) * (n - k + 1) FUNCTION = teller / noemer END FUNCTION SUB Oxalis (BYVAL hoek AS SINGLE, G AS Complex) EXPORT ' returns the normalized values for the function (-1 to + 1) ' for hoek = 0 the function start with real=1, imag =0 LOCAL Pol AS polar pol.ang = hoek pol.mag = COS(2!* hoek) P2C Pol, G END SUB FUNCTION Lemniskaat (BYVAL hoek AS SINGLE) EXPORT AS SINGLE ' the function normally is written: r^2 = a^2 . cos(2.hoek) ' if we take a = SQR(2) , the function becomes normalized. ' starts on 1 for hoek = 0 ' this functions changes hoek, such that it goes over and back in time! STATIC forward AS BYTE LOCAL Tmp AS SINGLE IF hoek > 6.283185307 THEN DO hoek -= 6.283185307 LOOP UNTIL hoek < 6.283185307 END IF IF hoek >= 3.141592654 THEN hoek = 6.283185307 - hoek forward = %False ELSE forward = %True END IF Tmp = SIN(2!*hoek) 'changed such that it starts on 0. COS(2! * hoek) IF Tmp < 0! THEN IF ISTRUE forward THEN FUNCTION = - SQR(ABS(Tmp)) ELSE FUNCTION = SQR(ABS(Tmp)) END IF ELSEIF Tmp > 0! THEN IF ISTRUE forward THEN FUNCTION = SQR(Tmp) ELSE FUNCTION = - SQR(Tmp) END IF ELSE FUNCTION = 0! END IF END FUNCTION ' averaging functions: FUNCTION AvgAri (BYVAL a!, BYVAL b!) EXPORT AS SINGLE FUNCTION = (a! + b!) / 2! END FUNCTION FUNCTION AvgGeo (BYVAL a!, BYVAL b!) EXPORT AS SINGLE LOCAL tmp AS SINGLE tmp = a! * b! IF tmp < 0 THEN FUNCTION = - SQR(ABS(tmp)) ELSE FUNCTION = SQR(tmp) END IF END FUNCTION FUNCTION AvgHar (BYVAL a!, BYVAL b!) EXPORT AS SINGLE IF a! + b! <> %False THEN FUNCTION = (2! * a! * b!)/ (a! + b!) ELSE FUNCTION = SGN(a!*b!) * 1E27 ' overflow !!! END IF END FUNCTION FUNCTION AvgMed (BYVAL a!, BYVAL b!) EXPORT AS SINGLE FUNCTION = SQR(((a!*a!) + (b!*b!))/ 2) ' was bug, we forgot the SQR ! corrected 20.02.2003 END FUNCTION ' averaging functions written for data acquisition buffers: ' parameters passed are : the data array itself and the number of elements over which to average ' the averaging starts at Ubound(array) - n FUNCTION AvgAriArr (dta() AS INTEGER, BYVAL n AS LONG) EXPORT AS SINGLE ' rekenkundig gemiddelde LOCAL i AS LONG LOCAL s AS SINGLE LOCAL sz AS LONG sz = UBOUND(dta) IF ISFALSE sz THEN FUNCTION = dta(0) : EXIT FUNCTION IF ISFALSE n THEN EXIT FUNCTION IF n = 1 THEN FUNCTION = dta(sz) : EXIT FUNCTION n = sz - n + 1 IF n < LBOUND(dta) THEN n = LBOUND(dta) FOR i = n TO sz s += dta(i) NEXT i FUNCTION = s / (sz-n+1) END FUNCTION FUNCTION AvgGeoArr (dta() AS INTEGER, BYVAL n AS LONG) EXPORT AS SINGLE ' meetkundig gemiddelde LOCAL i AS LONG LOCAL s AS DOUBLE LOCAL sz AS LONG sz = UBOUND(dta) IF ISFALSE sz THEN FUNCTION = dta(0) : EXIT FUNCTION n = sz - n + 1 IF n < LBOUND(dta) THEN n = LBOUND(dta) FOR i = n TO sz s *= dta(i) NEXT i FUNCTION = s ^ (1/(sz-n+1)) END FUNCTION FUNCTION AvgHarArr (dta() AS INTEGER, BYVAL n AS LONG) EXPORT AS SINGLE ' harmonisch gemiddelde ' do not call this function if the data can contain zeroes LOCAL i AS LONG LOCAL s AS DOUBLE LOCAL sz AS LONG sz = UBOUND(dta) IF ISFALSE sz THEN FUNCTION = dta(0) : EXIT FUNCTION n = sz - n + 1 IF n < LBOUND(dta) THEN n = LBOUND(dta) FOR i = n TO sz IF dta(i) <> %False THEN s += (1!/dta(i)) ELSE s += 1E27 ' overflow !!! Warning "Illegal value in data array in " + FUNCNAME$, 10000 END IF NEXT i IF s <> %False THEN FUNCTION = (sz-n+1) / s ELSE Warning "Overflow in g_indep.dll in " + FUNCNAME$, 10000 FUNCTION = %False END IF END FUNCTION FUNCTION AvgMedArr (dta() AS INTEGER, BYVAL n AS LONG) EXPORT AS SINGLE ' middelbare waarde ' calculates the mean of squares LOCAL i AS LONG LOCAL s AS SINGLE LOCAL sz AS LONG sz = UBOUND(dta) IF ISFALSE sz THEN FUNCTION = dta(0) : EXIT FUNCTION n = sz - n + 1 IF n < LBOUND(dta) THEN n = LBOUND(dta) FOR i = n TO sz s += (dta(i)^2) NEXT i FUNCTION = SQR(s / (sz-n+1) ) ' sqr added 14.02.2003, bug removed. END FUNCTION FUNCTION AvgSurf (dta() AS INTEGER, BYVAL n AS LONG) EXPORT AS SINGLE ' middelbare waarde ' voor oppervlakte bepaling in radarsystemen. LOCAL i AS LONG LOCAL s AS SINGLE LOCAL sz AS LONG sz = UBOUND(dta) IF ISFALSE sz THEN FUNCTION = dta(0) : EXIT FUNCTION n = sz - n + 1 IF n < LBOUND(dta) THEN n = LBOUND(dta) FOR i = n TO sz s += (dta(i)^2) NEXT i FUNCTION = s / (sz-n+1) ' here we do not square root, since the signal is proportional to the square root of ' the surface. So AvgSurf = AvgMedArr ^ 2 ' Thus the range of this function will be twice the range of the input. END FUNCTION ' statistiek FUNCTION AvgStatArr (dta() AS INTEGER, s AS SINGLE, BYVAL n AS LONG) EXPORT AS SINGLE LOCAL l AS SINGLE LOCAL i AS LONG LOCAL sz AS LONG LOCAL d AS SINGLE s = %False IF n < 2 THEN EXIT FUNCTION l = AvgAriArr (dta(),n) ' gewoon gemiddelde sz = UBOUND(dta) IF sz < 1 THEN FUNCTION = dta(0) : EXIT FUNCTION ' foutberekening: n = sz - n + 1 ' aantal metingen IF n < LBOUND(dta) THEN n = LBOUND(dta) FOR i = n TO sz d +=((dta(i) - l)^ 2) ' metode van de kleinste kwadraten. NEXT i 'changed 10.04.2010 to: s = SQR(d/n-1) ' d/n wanneer de populatie volledig is ' d/n-1 voor een random genomen populatie 'was: ' d = SQR(d/n) ' middelbare fout per meting '' s = d / SQR(n+1) ' standaard deviatie, middelbare fout in het gemiddelde ' statistics math memo: ' 68% of all values in dta are between l - s and l + s ' 95% are between l -2s and l +2s ' 97% are between l -3s and l +3s FUNCTION = l END FUNCTION ' library of functions for complex math. ' extensively used in electronics as well as ' in I/Q radar systems, where the output is a complex number in cartesian format. SUB AddCom (G1 AS Complex, G2 AS Complex, Gr AS Complex) EXPORT Gr.real = G1.real + G2.real Gr.imag = G1.imag + G2.imag END SUB SUB AddPol (G1 AS Polar, G2 AS Polar, Gr AS Polar) EXPORT ' eerst omzetten naar complex cartesisch DIM C1 AS Complex DIM C2 AS Complex DIM Cr AS Complex P2C G1, C1 P2C G2, C2 AddCom C1, C2, Cr C2P Cr, Gr END SUB FUNCTION AngCom# (C AS Complex) EXPORT LOCAL ang# ' this functions returns the phase angle in Rad for a complex number IF C.real <> 0 THEN ang# = ATN(C.imag / C.real) ELSE ang# = SGN(C.imag) * .1570796327 '3.141592654/ 2# END IF ' correction for quadrant-sign: 0-180 gr & 0-> -180gr IF C.real < 0 THEN ang# = 3.141592654 + ang# IF ang# > 3.141592654 THEN ang# = - 6.283185307 + ang# FUNCTION = ang# END FUNCTION SUB C2P (C AS Complex, P AS Polar) EXPORT ' converts cartesian complex numbers to polar coordinates P.mag = MagCom#(C) P.ang = AngCom#(C) END SUB SUB ConCom (G1 AS Complex, G2 AS Complex) EXPORT ' returns the conjugate (komplex toegevoegd getal) number G2.real = G1.real G2.imag = -G1.imag END SUB SUB ConPol (P1 AS Polar, P2 AS Polar) EXPORT ' complex conjugate in polar form DIM G1 AS Complex, G2 AS Complex P2C P1, G1 ConCom G1, G2 C2P G2, P2 END SUB FUNCTION Deg2Rad# (BYVAL degrees#) EXPORT Deg2Rad# = degrees# * 3.141592654 / 180# END FUNCTION SUB DivCom (G1 AS Complex, G2 AS Complex, Gr AS Complex) EXPORT ' eerst omzetten naar polair: DIM P1 AS Polar, P2 AS Polar, Pr AS Polar C2P G1, P1 C2P G2, P2 DivPol P1, P2, Pr P2C Pr, Gr ' ofwel, teller en noemer vermenigvuldigen met de conjugate: ' DIM ConjugateG2 AS Complex ' ConjugateG2.real = G2.real ' ConjugateG2.imag = -G2.imag ' DIM Teller AS Complex ' MulCom G1, ConjugateG2, Teller ' DIM Noemer AS Complex ' MulCom G2, ConjugateG2, Noemer END SUB SUB DivPol (G1 AS Polar, G2 AS Polar, Gr AS Polar) EXPORT Gr.mag = G1.mag / G2.mag Gr.ang = G1.ang - G2.ang END SUB SUB InvCom (G1 AS Complex, Gr AS Complex) EXPORT ' calculates 1/G1 DIM P1 AS Polar DIM P2 AS Polar P2.mag = 1# P2.ang = 0 DIM Pr AS Polar C2P G1, P1 DivPol P2, P1, Pr P2C Pr, Gr END SUB SUB InvPol (G1 AS Polar, Gr AS Polar) EXPORT ' calculates 1/G1 IF G1.mag <> 0 THEN Gr.mag = 1# / G1.mag ELSE Gr.mag = 1D+66 END IF Gr.ang = -G1.ang END SUB FUNCTION MagCom# (C AS Complex) EXPORT ' returns the magnitude or modulus of a complex number. ' This is the first part of the same number in polar form 'MagCom# = SQR((C.real ^ 2#) + (C.imag ^ 2#)) ' somewhat faster: MagCom# = SQR((C.real * C.real) + (C.imag * C.imag)) END FUNCTION SUB MulCom (G1 AS Complex, G2 AS Complex, Gr AS Complex) EXPORT ' eerst omzetten naar polair: DIM P1 AS Polar DIM P2 AS Polar DIM Pr AS Polar C2P G1, P1 C2P G2, P2 MulPol P1, P2, Pr P2C Pr, Gr END SUB SUB MulPol (G1 AS Polar, G2 AS Polar, Gr AS Polar) EXPORT Gr.mag = G1.mag * G2.mag Gr.ang = G1.ang + G2.ang END SUB FUNCTION OmegF# (BYVAL f#) EXPORT ' cirkelfrekwentie OmegF# = 6.283185307 * f# END FUNCTION SUB P2C (P AS Polar, C AS Complex) EXPORT C.real = P.mag * COS(P.ang) C.imag = P.mag * SIN(P.ang) END SUB SUB ParCom (G1 AS Complex, G2 AS Complex, G3 AS Complex) EXPORT ' 1/Z = 1/Z1 + 1/Z2 ==> Z = (Z1*Z2)/(Z1+Z2) DIM PsumC AS Complex, PsumP AS Polar, Pprod AS Polar, G1pol AS Polar, G2pol AS Polar, G3pol AS Polar AddCom G1, G2, PsumC C2P PsumC, PsumP: ' convert to polar C2P G1, G1pol C2P G2, G2pol MulPol G1pol, G2pol, Pprod DivPol Pprod, PsumP, G3pol P2C G3pol, G3 END SUB SUB ParPol (G1 AS Polar, G2 AS Polar, G3 AS Polar) EXPORT ' 1/Z = 1/Z1 + 1/Z2 ==> Z = (Z1*Z2)/(Z1+Z2) DIM Psum AS Polar, Pprod AS Polar AddPol G1, G2, Psum MulPol G1, G2, Pprod DivPol Pprod, Psum, G3 END SUB FUNCTION Rad2Deg# (BYVAL rad#) EXPORT Rad2Deg# = rad# * 57.29577951 ' rad# * 180# / Pi END FUNCTION SUB SubCom (G1 AS Complex, G2 AS Complex, Gr AS Complex) EXPORT Gr.real = G1.real - G2.real Gr.imag = G1.imag - G2.imag END SUB SUB SubPol (G1 AS Polar, G2 AS Polar, Gr AS Polar) EXPORT ' omzetten naar complex cartesisch DIM C1 AS Complex, C2 AS Complex, Cr AS Complex P2C G1, C1 P2C G2, C2 SubCom C1, C2, Cr C2P Cr, Gr END SUB SUB XcCar (BYVAL C#, BYVAL f#, Gr AS Complex) EXPORT Gr.real = 0 IF C# * f# > 0 THEN Gr.imag = -1# / (OmegF#(f#) * C#) ELSE Gr.imag = - 1D66 'Infinite END IF END SUB SUB XcPol (BYVAL C#, BYVAL f#, Gr AS Polar) EXPORT IF f# * C# > 0 THEN Gr.mag = 1# / (OmegF#(f#) * C#) ELSE Gr.mag = 1D66 ' Infinite END IF Gr.ang = - HalfPi '1.570796327 '(Pi /2!) END SUB SUB XlCar (BYVAL L#, BYVAL f#, Gr AS Complex) EXPORT Gr.real = 0 Gr.imag = OmegF#(f#) * L# END SUB SUB XlPol (BYVAL L#, BYVAL f#, Gr AS Polar) EXPORT Gr.mag = OmegF#(f#) * L# Gr.ang = HalfPi '1.570796327 ' 3.141592654/ 2! END SUB SUB DrawOxalis (BYVAL h AS LONG, p AS POINTL) EXPORT ' written to check our oxalis function in the Math library ' p hold the cooordinates of the start position to draw the graph LOCAL ppoint AS POINTL ' long LOCAL hDC AS LONG LOCAL G AS complex LOCAL i AS SINGLE hDC = GetDC(h) ' point.x = 210 ' point.y = 210 ScreenToClient h, p 'point DPtoLP hDC, p, 1 MoveToEx hDC, INT(p.x),INT(p.y), p FOR i = 0 TO 6.283185307 STEP 0.01 ' 0.062831853 Oxalis i, G 'point.x = INT(10 + 100 + (G.real * 100)) 'point.y = INT(210 + (G.imag * 100)) ppoint.x = INT(p.x + (G.real * 100)) ppoint.y = INT(p.y + (G.imag * 100)) ScreenToClient h, ppoint DPtoLP hDC, ppoint, 1 LineTo hDC, ppoint.x,ppoint.y NEXT i ReleaseDC h, hDC END SUB SUB DrawLemniskaat (BYVAL h AS LONG, p AS POINTL) EXPORT ' written to check our Lemniskaat function in the Math library ' p hold the coordinates of the start position to draw the graph LOCAL ppoint AS POINTL ' long LOCAL hDC AS LONG LOCAL G AS complex LOCAL i AS SINGLE LOCAL j AS SINGLE LOCAL t AS SINGLE hDC = GetDC(h) 'point.x = 10 'point.y = 210 ScreenToClient h, p DPtoLP hDC, p, 1 MoveToEx hDC, INT(p.x),INT(p.y), p FOR i = 0 TO 6.283185307 STEP 0.01 t = i j = Lemniskaat(t) ' note that this function changes t!!! ppoint.y = INT(p.y + (j * 100)) ppoint.x = INT(p.x + (t * 100)) ScreenToClient h, ppoint DPtoLP hDC, ppoint, 1 LineTo hDC, ppoint.x,ppoint.y NEXT i ReleaseDC h, hDC END SUB SUB DrawBarChart (BYVAL hw AS LONG, area AS rectl, arr() AS SINGLE) EXPORT ' use for statistic analysis. ' draw bar chart in window hw 'chart contains as much bars as elements in arr(), drawn in area 'values in arr() should be scaled [0, 1] , normalized. 'use hw = 0 for cleanup = release brushes 'note: we redraw entire window each time - remembering old values would make it impossible to have several ' barchart windows simultaneously... STATIC h AS INTEGER STATIC brush() AS LONG STATIC penblack AS LONG, pengrey AS LONG LOCAL bw%, Sp% LOCAL il AS LONG, H1% LOCAL hDC AS LONG LOCAL versize% LOCAL t AS DWORD LOCAL y AS LONG LOCAL nrbars AS LONG LOCAL hbt AS LONG LOCAL hpt AS LONG IF ISFALSE h THEN DIM brush(8) 'if we create & destroy brushes all the time windows fails creating new ones after a while.. h = 1 brush(0) = CreateSolidBrush(RGB(255, 64, 64)) brush(1) = CreateSolidBrush(RGB(255,64,127)) brush(2) = CreateSolidBrush(RGB(64,64,255)) brush(3) = CreateSolidBrush(RGB(64,127, 255)) brush(4) = CreateSolidBrush(RGB(64,255,64)) brush(5) = CreateSolidBrush(RGB(127,255,64)) brush(6) = CreateSolidBrush(RGB(255,255,127)) brush(7) = CreateSolidBrush(RGB(255,127,64)) brush(8) = CreateSolidBrush(RGB(255,127,255)) penblack = CreatePen(%PS_SOLID, 1, 0) pengrey = CreatePen(%PS_SOLID, 1, &H888888) END IF IF ISFALSE hw THEN FOR il = 0 TO 8 DeleteObject brush(il) h% = %false NEXT DeleteObject penblack DeleteObject pengrey EXIT SUB END IF nrbars = UBOUND(arr) - LBOUND(arr)+ 1 hDC = GetDC (hw) bw% = (area.nright - area.nleft)/nrbars 'x / nrbars sp% = bw% * 7/12 bw% = bw% '* 7/12 ' blank area: PatBlt hDC, Area.nleft,Area.nTop,Area.nRight, Area.nBottom ,%WHITENESS y = Area.nBottom - Area.nTop hbt = SelectObject(hDc, Brush(0)) FOR il = LBOUND(arr) TO UBOUND(arr) SelectObject hDC, Brush(il MOD 9) H1% = Area.nleft + (il * bw%) versize% = (Area.nBottom - (arr(il)* y)) Rectangle hDC, H1, Area.nBottom , H1%+Sp%, versize% NEXT il hpt = SelectObject (hDc, pengrey) FOR il = 1 TO 10 MoveToEx hDc, Area.nLeft, Area.nTop + il * y/10, BYVAL 0: LineTo hDc, Area.nRight, Area.nTop + il * y/10 NEXT SelectObject hDc, penblack MoveToEx hDc, Area.nLeft, Area.nTop + y/2, BYVAL 0: LineTo Hdc, Area.nRight, Area.nTop + y/2 SelectObject hDc, hpt Selectobject hDc, hbt ReleaseDC hw, hDC END SUB FUNCTION Lim2Pow2 (BYVAL value AS DWORD) EXPORT AS DWORD ' this function returns the largest number that is a power of 2 and that is smaller than the value passed. ' The function is used in combination with DFT's that ought to be calculated on numbers of samples that ' must be powers of 2. IF ISFALSE value THEN FUNCTION = %False : EXIT FUNCTION LOCAL n AS WORD n = LOG2(value) FUNCTION = EXP2(n) END FUNCTION FUNCTION NrSamp2Samptim (BYVAL NrSamp AS BYTE, BYVAL Sr AS SINGLE) EXPORT AS WORD ' returns the timeframe covered by a given number of samples in milliseconds. The sampling rate must be passed ' in Sr and is expressed in Hz IF ISFALSE Sr THEN FUNCTION = %False: EXIT FUNCTION FUNCTION = (NrSamp * 1000!) / Sr END FUNCTION FUNCTION SampTim2NrSamp (BYVAL SampTim AS WORD, BYVAL Sr AS SINGLE) EXPORT AS BYTE ' returns the number of samples required to cover a given timeinterval expressed in milliseconds. ' The samplingrate must be passed in Sr and is expressed in Hz IF ISFALSE Sr THEN FUNCTION = %False: EXIT FUNCTION FUNCTION = INT((Sr * Samptim) / 1000!) END FUNCTION FUNCTION Sinc (BYVAL angle AS SINGLE) EXPORT AS SINGLE IF angle <> 0 THEN FUNCTION = SIN(angle) / angle ELSE FUNCTION = %True END IF END FUNCTION FUNCTION Cosc (BYVAL angle AS SINGLE) EXPORT AS SINGLE LOCAL retval! IF angle <> 0 THEN retval! = COS(angle) / angle IF retval! > 1 THEN retval! = 1 IF retval! < -1 THEN retval! = -1 ELSE retval! = %False END IF FUNCTION = retval! END FUNCTION FUNCTION dBV2val (BYVAL dBv AS SINGLE) EXPORT AS WORD '[19.10.99] ' output range: 0 to + 2^16 ' dB range 0-90 dB returns 0-2^15. Note that sample data is bipolar! FUNCTION = (EXP10(ABS(dbv)/20)) / 2 END FUNCTION FUNCTION val2dBV (BYVAL value AS INTEGER) EXPORT AS SINGLE ' for 16 bit wave data, the range is 0-90dB IF ABS(value) >= 2 THEN FUNCTION = 20 * LOG10(ABS(value/2)) ELSE FUNCTION = %False END IF END FUNCTION FUNCTION ChargeCap (BYVAL RC AS SINGLE, BYVAL t AS SINGLE) EXPORT AS SINGLE ' for normalized 1 volt input. ' returns the voltage over a cap charged over R after time t ' if t=0, the result is 0 FUNCTION = 1! - EXP(1)^(-t/RC) END FUNCTION FUNCTION DisChargeCap (BYVAL RC AS SINGLE, BYVAL t AS SINGLE) EXPORT AS SINGLE ' if t=0 , the result is 1 FUNCTION = EXP(1) ^ (-t/RC) END FUNCTION FUNCTION RC_Cal (BYVAL prop AS SINGLE, BYVAL RC AS SINGLE) EXPORT AS SINGLE ' prop must be positive and normalised to 0 - 1 ' if prop = 0 , the result must be 'infinite' ' if prop = 1 , the result must be 0. ' this function returns the time it takes for the RC circuit to reach prop of its's initial charge. FUNCTION = - LOG(1-prop) * RC ' base e logarithm END FUNCTION ' wave generators: FUNCTION Triangle (BYVAL hoek AS SINGLE) EXPORT AS SINGLE ' periodic triangle wave generator function ' generates bipolar value -1 to +1 amplitude. Phase as for normal Sin() function IF hoek > 6.283185307 THEN DO hoek = hoek - 6.283185307 LOOP UNTIL hoek < 6.283185307 END IF IF hoek < 0 THEN DO hoek = hoek + 6.283185307 LOOP UNTIL hoek >= 0 END IF SELECT CASE hoek CASE 0 TO HalfPi '1.570796327 'Pi/2 3.141592654 FUNCTION = hoek / 1.570796327 CASE HalfPi TO 4.71238898 '3/2 Pi FUNCTION = 1! - ((hoek - HalfPi)/HalfPi) CASE 4.71238898 TO Pi2 '6.283185307 FUNCTION = -1! + ((hoek - 4.71238898)/HalfPi) END SELECT END FUNCTION FUNCTION Saw (BYVAL hoek AS SINGLE, BYVAL DC AS SINGLE) EXPORT AS SINGLE ' 02.12.1999 ' this function starts at -1 for angle 0 input. LOCAL maxval AS SINGLE LOCAL sw AS SINGLE ' limit angle to a single full period: IF hoek > Pi2 THEN DO hoek -= Pi2 '6.283185307 LOOP UNTIL hoek < Pi2 '6.283185307 END IF IF hoek < 0 THEN DO hoek += Pi2 '6.283185307 LOOP UNTIL hoek >= 0 END IF maxval = Pi2 * DC '6.283185307 * DC SELECT CASE hoek CASE < maxval sw = hoek / maxval ' 0-1 CASE >= maxval sw = 1! - ((hoek - maxval) / maxval) ' 1-0 END SELECT ' convert to bipolar: FUNCTION = (sw + sw) - 1! ' avoid multiply: (sw*2)-1 END FUNCTION FUNCTION Square (BYVAL hoek AS SINGLE, BYVAL DC AS SINGLE) EXPORT AS INTEGER ' 02.12.1999 ' squarewave generator ' this function starts at 1 for angle 0 input. ' it returns either -1 or +1. LOCAL maxval AS SINGLE ' limit angle to a single full period: IF hoek > 6.283185307 THEN DO hoek = hoek - 6.283185307 LOOP UNTIL hoek < 6.283185307 END IF IF hoek < 0 THEN DO hoek = hoek + 6.283185307 LOOP UNTIL hoek >= 0 END IF maxval = 6.283185307 * DC SELECT CASE hoek CASE < maxval FUNCTION = %True CASE ELSE FUNCTION = %NotFalse END SELECT END FUNCTION FUNCTION ReduceAngle (BYVAL hoek AS SINGLE, BYVAL param AS BYTE) EXPORT AS SINGLE LOCAL Revolutions AS LONG Revolutions = FIX(hoek / Pi2) SELECT CASE param CASE %Unipolar ' return hoek as a positive angle between 0 and 2*Pi FUNCTION = hoek - (Revolutions * Pi2) CASE %Bipolar ' returns hoek as a signed value between 0 and Pi (positive) , Pi and 2*Pi=0 negative IF hoek - (Revolutions * Pi2) < Pi THEN FUNCTION = hoek - (Revolutions * Pi2) ELSE FUNCTION = (hoek - ((Revolutions * Pi2) + Pi)) END IF END SELECT END FUNCTION FUNCTION Gray (BYVAL n AS DWORD) EXPORT AS DWORD LOCAL b AS DWORD b = n SHIFT RIGHT b,1 FUNCTION = n XOR b END FUNCTION FUNCTION Kummer (BYVAL x AS SINGLE, BYVAL y AS SINGLE, BYVAL z AS SINGLE, BYVAL mu AS SINGLE) EXPORT AS SINGLE ' mu^2 should be choosen away from 1/3, 1 and 3. ' with mu = 13/10 you get the picture as on the cover of Linux 6.1 ' ref.: http://www-groups.dcs.st-and.ac.uk/~history/mathematicians/kummer.html LOCAL lambda AS SINGLE LOCAL p AS SINGLE LOCAL q AS SINGLE LOCAL r AS SINGLE LOCAL s AS SINGLE LOCAL noemer AS SINGLE STATIC ctrlwrd% noemer = 3-(mu^2) IF noemer <> 0! THEN lambda = ((3 * (mu^2)) - 1 ) / noemer ELSE Warning "Illegal value for mu-parameter passed in Kummer function in " + FUNCNAME$, 10000 ctrlwrd% = &B0001001100111111 ! fldcw ctrlwrd% ; guarantee preservation of extended precision - windows bug. FUNCTION = %False EXIT FUNCTION END IF p = 1 - z - (SQR(2) * x) q = 1 - z + (SQR(2) * x) r = 1 + z - (SQR(2) * y) s = 1 + z + (SQR(2) * y) FUNCTION = (((x^2) + (y^2) + (z^2) + (mu^2))^2) - (lambda * p * q * r * s) END FUNCTION FUNCTION Str2Fact (factorial AS STRING) EXPORT AS QUAD ' the input string had a single byte for each exponent. These bytes have hex format. ' Thus, the exponents are limited to the range 0-F, or, 0-15 ' used in Wave2Enveloppe. LOCAL j AS QUAD LOCAL i AS DWORD LOCAL cnt AS WORD STATIC tmp AS BYTE STATIC Primes() AS LONG STATIC tog AS BYTE LOCAL ctrlwrd% IF ISFALSE tog THEN DIM Primes(12) AS STATIC LONG Primes(1) = 1 ' not required. Primes(2) = 2 Primes(3) = 3 Primes(4) = 5 Primes(5) = 7 Primes(6) = 11 Primes(7) = 13 Primes(8) = 17 Primes(9) = 19 Primes(10) = 23 Primes(11) = 29 Primes(12) = 31 tog = %True END IF IF LEN(factorial) > 12 THEN 'Warning "Overflow in Str2Fact" & CHR$(12,13), 2000 Warning "Overflow - g_indep in " + FUNCNAME$, 10000 ctrlwrd% = &B0001001100111111 ! fldcw ctrlwrd% ; guarantee preservation of extended precision - windows bug. FUNCTION = -1 EXIT FUNCTION END IF j = 1 cnt = 1 FOR i = 1 TO LEN(factorial) cnt = 1 + LEN(factorial) - i ' i = 1 ===> cnt = lsb tmp = VAL("&H" & (MID$(factorial,cnt,1))) j = j * INT(Primes(i) ^ tmp) NEXT i FUNCTION = j END FUNCTION ' transform functions: SUB HanningWindow (Shape!()) EXPORT ' Version for single precision. ' to avoid repeated calls to this function, you should whenever you need ' the hanning window applied to your data, call it once with an array ' prepared as MAT arr!() = CON , this sets all elements to 1, and passed ' to this procedure it will return a normalized hanning array. In your code ' you can than apply the hanning window ' simply by MAT MyDataArray() = MyDataArray() * Arr!() LOCAL siz AS LONG LOCAL HalfSize AS LONG LOCAL QuartSize AS LONG LOCAL i AS LONG LOCAL hoekincrement!, hoek!, faktor1!, faktor2! siz = UBOUND(Shape!) hoekincrement! = Pi2 / siz hoek! = Pi HalfSize = siz / 2 QuartSize = siz / 4 DO faktor1! = (COS(hoek!) + 1) / 2 faktor2! = 1 - faktor1! Shape!(i) *= faktor1! Shape!(siz - i) *= faktor1! Shape!(HalfSize - i) *= faktor2! Shape!(HalfSize + i) *= faktor2! INCR i hoek! += hoekincrement! LOOP UNTIL i = QuartSize ' now we missed only Shape!(QuartSize) faktor1! = (COS(hoek!) + 1) / 2 Shape!(QuartSize) *= faktor1! END SUB SUB HanningWindow_Dbl (Shape#()) EXPORT ' would be much faster using a MAT Arr() function if we calculate the ' factors once beforehand. LOCAL siz AS LONG LOCAL HalfSize AS LONG LOCAL QuartSize AS LONG LOCAL i AS LONG LOCAL hoekincrement, hoek, faktor1, faktor2 AS DOUBLE siz = UBOUND(Shape#) hoekincrement = Pi2 / siz hoek = Pi HalfSize = siz / 2 QuartSize = siz / 4 DO faktor1 = (COS(hoek) + 1) / 2 faktor2 = 1 - faktor1 Shape#(i) *= faktor1 Shape#(siz - i) *= faktor1 Shape#(HalfSize - i) *= faktor2 Shape#(HalfSize + i) *= faktor2 INCR i hoek += hoekincrement LOOP UNTIL i = QuartSize ' now we missed only Shape!(QuartSize) faktor1 = (COS(hoek) + 1) / 2 Shape#(QuartSize) *= faktor1 END SUB SUB HanningWindow_dbl2sng (InArray#(), OutArray!()) EXPORT ' as the prevouw procedure, Hanningwindow modifies the array passed, ' we made this alternative, returning a new array of singles. ' input array and output array must have the same size! LOCAL siz AS LONG LOCAL HalfSize AS LONG LOCAL QuartSize AS LONG LOCAL i AS LONG LOCAL hoekincrement!, hoek!, faktor1!, faktor2! siz = UBOUND(InArray#) hoekincrement! = Pi2 / siz hoek! = Pi HalfSize = siz / 2 QuartSize = siz / 4 DO faktor1! = (COS(hoek!) + 1) / 2 faktor2! = 1 - faktor1! OutArray!(i) = InArray#(i) * faktor1! OutArray!(siz - i) = InArray#(siz - i) * faktor1! OutArray!(HalfSize - i) = InArray#(HalfSize - i) * faktor2! OutArray!(HalfSize + i) = InArray#(HalfSize + i) * faktor2! INCR i hoek! = hoek! + hoekincrement! LOOP UNTIL i = QuartSize ' now we missed only InArray!(QuartSize) faktor1! = (COS(hoek!) + 1) / 2 OutArray!(QuartSize) = InArray#(QuartSize) * faktor1! END SUB SUB HanningWaveWindow (BYVAL lwi AS INTEGER PTR, BYVAL s AS DWORD) EXPORT ' local ctrlwrd% ' this procedure returns the wave array passed with a hanning window applied to the data... Warning "Hanningwavewindow: In the works...", 5000 ',FUNCNAME$ ' ctrlwrd% = &B0001001100111111 ' ! fldcw ctrlwrd% ; guarantee preservation of extended precision - windows bug. END SUB SUB TriangleWindow (Shape!()) EXPORT LOCAL maat%, HalfSize% , i% LOCAL increment!, faktor! maat% = UBOUND(Shape!) HalfSize% = maat% / 2 increment! = 1! / HalfSize% faktor! = 0 FOR i% = 0 TO HalfSize% Shape!(i%) *= faktor! Shape!(maat% - i%) *= faktor! faktor! += increment! NEXT i% END SUB SUB BetaWindow (BYREF Sp!(), BYVAL a AS SINGLE, BYVAL b AS SINGLE) EXPORT 'we presume that the lbound of the array is 0 'returns a normalized Sp() array 'if you want a normalized array, with an empty array on input, use first: ' MAT Sp!()= CON to make all elements equal to 1 prior to calling this function ' The function was originally written to be used in conjuntion with the DFT procedures. LOCAL Siz AS LONG LOCAL i AS LONG LOCAL normval AS SINGLE LOCAL x AS SINGLE 'LOCAL ctrlwrd% normval = ((a^a) * (b^b))/((a + b)^(a + b)) IF ISFALSE normval THEN Warning "Parameters beyond range passed in module [g_indep] in " + FUNCNAME$, 10000 'ctrlwrd% = &B0001001100111111 '! fldcw ctrlwrd% ; guarantee preservation of extended precision - windows bug. EXIT SUB END IF Siz = UBOUND(Sp!) IF ISFALSE Siz THEN EXIT SUB FOR i = 0 TO Siz Sp!(i) *= GetBeta (a,b, i/siz) NEXT i END SUB FUNCTION GetBeta (BYVAL a AS SINGLE, BYVAL b AS SINGLE, BYVAL x AS SINGLE) EXPORT AS SINGLE STATIC oa AS SINGLE STATIC ob AS SINGLE STATIC normval AS SINGLE IF (a <> oa) OR (b <> ob) THEN normval = ((a^a) * (b^b))/((a + b)^(a + b)) oa = a ob = b IF ISFALSE normval THEN Warning "Parameters beyond range passed in [g_indep] in " + FUNCNAME$, 10000 FUNCTION = %False EXIT FUNCTION END IF END IF FUNCTION = ((x ^ a) * (1! - x) ^b)/normval END FUNCTION SUB DFT (Samp!(), Sp!()) EXPORT ' takes ca. 9ms on Miel on 256 samples. ' coding for single precision numbers ' improved for speed 04.05.2010 gwr LOCAL x AS LONG LOCAL i AS LONG LOCAL maat AS LONG LOCAL halfsize AS LONG LOCAL kfaktor!, kt!, Rl!, Im!, g!, pw! DIM a(UBOUND(Samp!)) AS LOCAL SINGLE maat = (UBOUND(Samp!)) + 1 ' must be a power of 2 halfsize = maat / 2 MAT a() = (1.0/maat) * Samp!() ' moved out of the loop for speedup REDIM Sp!(halfsize) kfaktor! = Pi2 / maat RESET i DO kt! = kfaktor! * i ' could become a lookup kt!(i), 0-halfsize RESET Rl!, Im! FOR x = 0 TO maat - 1 g! = kt! * x ' would become g! = kt(i) * x Rl! += (a(x) * COS(g!)) Im! -= (a(x) * SIN(g!)) NEXT x pw! = (Rl! * Rl!) + (Im! * Im!) ' sum of squares IF pw! THEN Sp!(i) = 2 * SQR(pw!) INCR i LOOP UNTIL i = halfsize ' the second half would be empty END SUB SUB DFT256_Dbl (Samp#(), Sp#()) EXPORT 'modified for speed 04.05.2010 - gwr '09.05.2010: optimized for use in Doppler gesture analysis context ' hanning window always applied. ' this version uses a lookup for the cos/sin factors (64kB) ' now works fine... ' coding can still be improved for compactness and speed ' In the interest of performance, no error checking is performed here. ' The arrays passed must be of the expected size and type!!! (255 and 127) LOCAL x, i AS LONG LOCAL Rl#, Im#, pw# STATIC tog AS DWORD DIM a(255) AS LOCAL DOUBLE DIM arr(255) AS LOCAL DOUBLE DIM cosg(255, 127) AS STATIC DOUBLE DIM sing(255, 127) AS STATIC DOUBLE DIM han(255) AS STATIC DOUBLE IF ISFALSE tog THEN ' prepare a hanning window data array of 256 elements: MAT arr()= CON(1) ' set all elements to 1 HanningWindow_Dbl arr() ' now apply the window to obtain a normalised array MAT han() = (3.90625E-3) * arr() FOR i = 0 TO 127 FOR x = 0 TO 255 cosg(x,i) = COS(0.0245436926 * i * x) '32kB sing(x,i) = SIN(0.0245436926 * i * x) '32kB NEXT x NEXT i RESET i 'we forgot this and had a big crash... tog = %True END IF MAT a() = Samp#() * Han() 'apply the hanning window to the sample array DO RESET Rl# , Im# FOR x = 0 TO 255 Rl# += (a(x) * cosg(x, i)) Im# -= (a(x) * sing(x, i)) NEXT x pw# = (Rl# * Rl#) + (Im# * Im#) ' sum of squares Sp#(i) = 2 * SQR(pw#) ' power spectrum INCR i LOOP UNTIL i = 128 END SUB SUB Cepstrum128_Dbl (Sp#(), Ce#()) EXPORT 'this procedure accepts a power spectrum with 128 bins as input 'and returns a cepstrum containing 64 bins in the quefrency domain as output 'the input normally comes from the DFT256_Dbl procedure 'no windowing is applied to the input spectrum. 'to be checked and debugged. Not sure whether this is mathematically correct 'we try using this to discover periodicity in the 4 second slow movement buffers (doppler namuda code) LOCAL x, i AS LONG LOCAL Rl#, Im#, pw# STATIC tog AS DWORD DIM a(127) AS LOCAL DOUBLE 'for lookup tables: DIM cosg(127, 63) AS STATIC DOUBLE DIM sing(127, 63) AS STATIC DOUBLE IF ISFALSE tog THEN FOR i = 0 TO 63 FOR x = 0 TO 127 cosg(x,i) = COS(0.04908738521 * i * x) ' 0.0490... = 2 Pi / 128 sing(x,i) = SIN(0.04908738521 * i * x) NEXT x NEXT i RESET i tog = %True END IF ' first take the log of the input spectrum: ' MAT a() = (1/128) * Sp#() ' can we skip this as it merely changes scale. FOR x = 0 TO 127 ' a(x) = ABS(LOG10(Sp#(x))) 'LOG alone gives all negative values in a(x) , does'nt seem to make a difference a(x) = LOG10(Sp#(x)) 'shifting by stating a(x) = 20 + LOG(Sp#(x)) does make no difference, although everything becomes positive 'potential error here: avoidance of LOG 0 !!! NEXT x 'ref. Curtis Roads, p.516 FOR i = 0 TO 63 RESET Rl# , Im# FOR x = 0 TO 127 Rl# += (a(x) * cosg(x, i)) Im# -= (a(x) * sing(x, i)) 'do we need the imag. part here? NEXT x 'for normal DFT: pw# = (Rl# * Rl#) + (Im# * Im#) ' sum of squares Ce#(i) = 2 * SQR(pw#) ' power spectrum 'taking only the real values: ' Ce#(i) = Rl# * Rl# ' Rl# alone gives a bipolar output NEXT i END SUB SUB Cepstrum128_Dbl_Inv (Sp#(), Ce#()) EXPORT 'this procedure accepts a power spectrum with 128 bins as input 'and returns a cepstrum containing 64 bins in the quefrency domain as output 'the input normally comes from the DFT256_Dbl procedure 'no windowing is applied to the input spectrum. 'to be checked and debugged. Not sure whether this is mathematically correct 'we try using this to discover periodicity in the 4 second slow movement buffers (doppler namuda code) 'here we apply an inverted DFT 'no good as yet!!! 'the resulting waveform is in principe twice the size of the spectrum, but it is symmetrical and thus 'consists of the mirrored halves. This fundamental frequency is the result of the windowing. 'here we seem to need only half of the quefrency convolution, and thus the size of Ce# should be the same 'as that of the input spectrum. (untrue: halfsize is enough) LOCAL x, i AS LONG LOCAL Rl# ', Im#, pw# STATIC tog AS DWORD DIM a(127) AS LOCAL DOUBLE 'for inverted DFT we have: DIM sing(127,63) AS STATIC DOUBLE IF ISFALSE tog THEN FOR x = 0 TO 63 FOR i = 0 TO 127 sing(i,x) = SIN(0.04908738521 * i * x) '2 Pi / 128 NEXT x NEXT i RESET i tog = %True END IF ' first take the log of the input spectrum: ' MAT a() = (1/128) * Sp#() ' can we skip this as it merely changes scale. FOR x = 0 TO 127 ' a(x) = ABS(LOG10(Sp#(x))) 'LOG alone gives all negative values in a(x) , does'nt seem to make a difference a(x) = LOG10(Sp#(x)) NEXT x 'ref. Curtis Roads, p.516 FOR x = 0 TO 63 'if we go to 127, we get a symmetrical wave RESET Rl# FOR i = 0 TO 127 'this always has to go through the complete spectrum Rl# += (a(i) * sing(i,x)) NEXT i Ce#(x) = Rl# * Rl# 'ABS(Rl#) 'quefrencies NEXT x END SUB SUB DFT_Dbl (Samp#(), Sp#()) EXPORT 'modified for speed 04.05.2010 - gwr LOCAL x AS LONG LOCAL i AS LONG LOCAL maat AS LONG LOCAL halfsize AS LONG LOCAL kfaktor#, kt#, Rl#, Im#, g#, pw# DIM a(UBOUND(Samp#)) AS LOCAL DOUBLE maat = (UBOUND(Samp#)) + 1 ' must be a power of 2 halfsize = maat / 2 'REDIM Sp#(halfsize) gives a problem with absolute arrays dimensioned with DIM AT and passed byref 'moreover: halfsize-1 is large enough, since we have no f-bin at Sp#(halfsize) kfaktor# = Pi2 / maat ' ATN(1)* 8# / maat MAT a() = (1.0/maat) * Samp#() 'for speedup placed outside the loop 04.05.2010 DO kt# = kfaktor# * i RESET Rl# , Im# FOR x = 0 TO maat - 1 g# = (kt# * x) Rl# += (a(x) * COS(g#)) Im# -= (a(x) * SIN(g#)) NEXT x pw# = (Rl# * Rl#) + (Im# * Im#) ' sum of squares IF pw# THEN Sp#(i) = 2 * SQR(pw#) INCR i LOOP UNTIL i = halfsize ' the second half will be empty END SUB SUB DFT_Long (Samp() AS LONG, Spektrum() AS LONG) EXPORT ' long integer version of above procedure. ' It's 10% faster, but does not require a conversion to float on entry. LOCAL i AS LONG LOCAL x AS LONG LOCAL Rl! LOCAL Im! LOCAL g! STATIC maat AS LONG STATIC halfsize AS LONG LOCAL pw! IF UBOUND(Samp) + 1 <> maat THEN ' calculate lookup only once... maat = UBOUND(Samp) + 1 halfsize = maat / 2 REDIM kt(halfsize) AS STATIC SINGLE FOR i = 0 TO halfsize kt(i) = (i * Pi2) / maat NEXT i 'REDIM spektrum(maat-1) REDIM Spektrum(halfsize) AS LONG REDIM a(maat-1) AS STATIC SINGLE END IF 'mat a() = (1/maat) * Samp() - kan niet omdat a() en samp() een verschillend type zijn. FOR i = 0 TO maat -1 a(i) = Samp(i) / maat NEXT i i = 0 DO Rl! = 0 Im! = 0 FOR x = 0 TO maat - 1 IF a(x) THEN g! = kt(i) * x Rl! += (a(x) * COS(g!)) Im! -= (a(x) * SIN(g!)) END IF NEXT x pw! = (Rl! * Rl!) + (Im! * Im!) ' sum of squares IF pw! THEN Spektrum(i) = 2 * SQR(pw!) INCR i LOOP UNTIL i = halfsize ' the second half will be empty END SUB SUB DFT_Int (Samp() AS INTEGER, Spektrum() AS INTEGER) EXPORT ' short integer version of above procedure. ' It's 10% faster, but does not require a conversion to float on entry. LOCAL i AS LONG LOCAL x AS LONG LOCAL Rl! LOCAL Im! LOCAL g! STATIC maat AS LONG STATIC halfsize AS LONG LOCAL pw! IF UBOUND(Samp) + 1 <> maat THEN ' calculate lookup only once... maat = UBOUND(Samp) + 1 halfsize = maat / 2 REDIM kt(halfsize) AS STATIC SINGLE FOR i = 0 TO halfsize kt(i) = (i * Pi2) / maat NEXT i 'REDIM spektrum(maat-1) REDIM spektrum(halfsize) AS INTEGER REDIM a(maat-1) AS STATIC SINGLE END IF 'MAT a() = (1.0/maat) * Samp() gaat niet want de arrays hebben een verschillend type... FOR i = 0 TO maat -1 a(i) = Samp(i) / maat NEXT i i = 0 DO Rl! = 0 Im! = 0 FOR x = 0 TO maat - 1 IF a(x) THEN g! = kt(i) * x Rl! += (a(x) * COS(g!)) Im! -= (a(x) * SIN(g!)) END IF NEXT x pw! = (Rl! * Rl!) + (Im! * Im!) ' sum of squares IF pw! THEN Spektrum(i) = 2 * SQR(pw!) INCR i LOOP UNTIL i = halfsize ' the second half will be empty END SUB SUB InvDFT (Sp!(), Smp!()) EXPORT LOCAL siz AS LONG LOCAL i AS LONG LOCAL x AS LONG LOCAL kfaktor!, kt!, Rl!, g!, a! ' transforms frequency domain into timedomain ' Time assumed to be always 1 second for the Samp!() length. siz = UBOUND(Sp!) + 1 ' REDIM Smp!(UBOUND(Sp!)) ' moet tweemaal zo groot zijn. REDIM Smp!((siz * 2) -1) kfaktor! = Pi2 / siz DO kt! = kfaktor! * i RESET Rl!, x DO g! = (kt! * x) a! = Sp!(x) IF a! > 0 THEN a! = SQR(a!) Rl! += (a! * SIN(g!)) ' only real values END IF INCR x LOOP UNTIL x >= siz / 2: ' half size is enough ' IF Rl! THEN Smp!(i) = Rl! ' END IF INCR i ' counts frequency bands LOOP UNTIL i > siz - 1 END SUB FUNCTION CreateDefaultIniFile (f AS STRING) EXPORT AS LONG ' the string constants used here are in gmt_kons.bi. ' This procedure should only be called if the application does not find a valid .INI file. LOCAL nr AS LONG LOCAL i% LOCAL App AS ApplicationType LOCAL mpName AS ASCIIZ * %MAX_COMPUTERNAME_LENGTH + 1 DIM Spec(1, 16) AS LOCAL SINGLE 'INTEGER - was bug before 25.12.2004 !!! DIM FuzInt(11) AS LOCAL SINGLE DIM FuzMel(-11 TO 11) AS LOCAL SINGLE ' DIM AudioFader(0 TO 1) AS AudioFaderType ' to be changed !!! App.komposduur = 3600 ' in seconds - one hour default App.tempo = 60 ' in M.M. App.globton = %Null App.DebugTaskNr = %False App.WriteSeqScoreTaskNr = 5 ' dll tasks. App.ReadSeqScoreTaskNr = 10 App.ShowGlobalHarmonyTaskNr = 14 App.GlobalHarmonyTaskNr = 15 App.Audiofilepath = "C:\audiofiles\" App.Midifilepath = "C:\midifiles\" App.Datafilepath = "C:\b\pb\gmt\" GetComputerName mpName,SIZEOF(mpName) ' WinApi function mpname = MCASE$(TRIM$(mpname)) IF f = "" THEN f = TRIM$(mpname) & "_" & TRIM$($INI) nr = FREEFILE OPEN f FOR OUTPUT AS #nr PRINT #nr," Default Initialisation file created by {g_indep.dll} module" PRINT #nr,"-----------------------------------------------------------------------" PRINT #nr,"This file was created for the PC named "; mpname PRINT #nr, SPACE$(10) PRINT #nr, "[WIN_DIR],", GetWindir PRINT #nr, SPACE$(10) PRINT #nr, $ApplicationData PRINT #nr, $ApDur & ",", App.komposduur '"KOMPOSITIEDUUR,";App.komposduur PRINT #nr, $ApGT & ",", App.globton '"TONAL_CENTER,";App.globton PRINT #nr, $ApTempo & ",", App.tempo '"TEMPO,";App.tempo PRINT #nr, $Pitch & ",", 440! '"DIAPASON,";440 PRINT #nr, "AUTEUR," , "-----" PRINT #nr, "PROJEKTNAAM," , "UnNamed" PRINT #nr, $ApAudiofilepath, App.Audiofilepath PRINT #nr, $ApMidifilepath, App.midifilepath PRINT #nr, $ApDataFilePath, App.datafilepath PRINT #nr, $ApplicationDataEnd PRINT #nr, SPACE$(10) PRINT #nr, $Cockpit PRINT #nr, $NrUD & ",", 18 '%False creates 1 ud controller PRINT #nr, $NrSL & ",", 2 ' such that the demo works... PRINT #nr, $NrTOG & ",",12 PRINT #nr, $NrBUT & ",",12 PRINT #nr, $Cockpit_Labels PRINT #nr, $Cockpit_End PRINT #nr, SPACE$(10) PRINT #nr, SPACE$(10) PRINT #nr, $TSK_START PRINT #nr, $TSK_WRISEQ; ",", App.WriteSeqScoreTaskNr '"WRITE_SEQUENCE_TASK_NUMBER," PRINT #nr, $TSK_RDSEQ; ",", App.ReadSeqScoreTaskNr '"READ_SEQUENCE_TASK_NUMBER, " PRINT #nr, $TSK_SHOWHAR; ",", App.ShowGlobalHarmonyTaskNr '"SHOW_GLOBAL_HARMONY_TASK_NUMBER," PRINT #nr, $TSK_GLOBHAR; ",", App.GlobalHarmonyTaskNr '"GLOBAL_HARMONY_TASK_NUMBER," PRINT #nr, $TSK_DEBUG; ",", App.DebugTaskNr '"DEBUG_TASK_NUMBER," PRINT #nr, $TSK_END PRINT #nr, SPACE$(10) ' should also be removed: ' PRINT #nr, $Audio_Start ' PRINT #nr, $CD_VOL_L;",", &HFFF0 ' PRINT #nr, $CD_VOL_R;",", &HFFF0 ' PRINT #nr, $AX_VOL_L;",", &HFFF0 ' PRINT #nr, $AX_VOL_R;",", &HFFF0 ' PRINT #nr, $Audio_End PRINT #nr, SPACE$(10) PRINT #nr, $Flags_Start PRINT #nr, $FlagScore; ", 0" PRINT #nr, $FlagMidi; ", 0" PRINT #nr, $FlagHarm; ", 0" PRINT #nr, $ApFlga; ", 1" ', 1 '"AUTOPATCH,"; 1 ' if set, tasks will send midi patch info when switched on. PRINT #nr, $FlagSysEx; ", 0" PRINT #nr, $FlagSysEx2File; ",0" PRINT #nr, $FlagSysEx2Array; ",0" PRINT #nr, $Flags_End ' The information sent, is in Task().channel, .patch, level, .pan... PRINT #nr, SPACE$(10) PRINT #nr, $Midi_Start PRINT #nr, "nrgroups, 1" PRINT #nr, "0, GM" PRINT #nr, "GROUP 0" PRINT #nr, "nrdevs, 1" PRINT #nr, "0, Generic, -1, 0, 0" PRINT #nr, $Midi_End ' optional entries for sensors and specific robots: PRINT #nr, "[CQT]" PRINT #nr, "[INPUT_PORT],0" PRINT #nr, "[OUTPUT_PORT],4" PRINT #nr, "[INPUT_CHANNEL], 0" PRINT #nr, "[OUTPUT_CHANNEL], 0" PRINT #nr, "[LOWTES], 32" PRINT #nr, "[HIGHTES], 93" PRINT #nr, "[CC33_ONSET], 75" PRINT #nr, "[CC34_OFFSET], 85" PRINT #nr, "[CC35_META_Q_C1_SENS], 64" PRINT #nr, "[CC36_META_Q_C2_SENS], 64" PRINT #nr, "[CC37_META_S_C5], 64" PRINT #nr, "[CC38_META_S_C7], 64" PRINT #nr, "[CC39_META_S_C9], 64" PRINT #nr, "[CC40_TRANS], 64" PRINT #nr, "[CC41_HIGHTES], 93" PRINT #nr, "[CC42_LOWTES], 32" PRINT #nr, "[CQT_END]" PRINT #nr, "' tuning for Hurdy: (19.10.2004 - La- Re stemming)" PRINT #nr, "[HURDY_TUNING]" PRINT #nr, "33" PRINT #nr, "50" PRINT #nr, "[HURDY_TUNING_END]" PRINT #nr, "[HY1]" PRINT #nr, "[INPUT_PORT], 0" PRINT #nr, "[INPUT_CHANNEL], 0" PRINT #nr, "[HY1_END] PRINT #nr, "[AXE3]" PRINT #nr, "[INPUT_PORT], 0" PRINT #nr, "[INPUT_CHANNEL], 0" PRINT #nr, "[AXE3_END] PRINT #nr, SPACE$(10) PRINT #nr, "/Data section used in the harmony library/" PRINT #nr, "/Fuzzy values for interval dissonance/" PRINT #nr, $FuzHar 'GetDefaultFuzzyHarValues FuzInt() ' in gmt_lib.dll FuzInt(0)= 0 FuzInt(1)= 1 FuzInt(2)= 0.75 FuzInt(3)= 0.4 FuzInt(4)= 0.3 FuzInt(5)= 0.05 FuzInt(6)= 0.85 FuzInt(7)= 0.05 FuzInt(8)= 0.27 FuzInt(9) = 0.32 FuzInt(10)= 0.67 FuzInt(11)= 0.91 FOR i% = 0 TO UBOUND(FuzInt) PRINT #nr, STR$(i%) & ", " & STR$(FuzInt(i%)) NEXT i% PRINT #nr, $FuzHar_End PRINT #nr, SPACE$(10) PRINT #nr, "/Fuzzy values for melodic property strength/" PRINT #nr, $FuzMel 'GetDefaultFuzzyMelValues FuzMel() FuzMel(0)= .02 FuzMel(1)= .06 FuzMel(-1)= .06 FuzMel(2)= .04 FuzMel(-2)= .04 FuzMel(3)= .1 FuzMel(-4)= .1 FuzMel(4)= .2 FuzMel(-3)= .2 FuzMel(5)= .35 FuzMel(-7)= .35 FuzMel(6)= .8 FuzMel(-6)= .8 FuzMel(7)= .08 FuzMel(-5)= .08 FuzMel(8)= .7 FuzMel(-9)= .7 FuzMel(9)= .6 FuzMel(-8)= .6 FuzMel(10)= .7 FuzMel(-10)= .7 FuzMel(-11)= .81 FuzMel(11)= .81 PRINT #nr, STR$(%False) & ", " & STR$(FuzMel(%False)) FOR i% = 1 TO UBOUND(FuzMel) PRINT #nr, STR$(-i%) & ", " & STR$(FuzMel(-i%)) PRINT #nr, STR$(i%) & ", " & STR$(FuzMel(i%)) NEXT i% PRINT #nr, $FuzMel_End PRINT #nr, SPACE$(10) PRINT #nr, "/Spectral data used in spectral harmony calculation/" PRINT #nr, $Spectrum i% = %False DO Spec(0, i%) = i% + 1! Spec(1, i%) = 1! / (i% + 1!) ' must be single - was bug before 25.12.2004 PRINT #nr, STR$(i%) & ", " & STR$(%False) & ", " & STR$(Spec(0,i%)) PRINT #nr, STR$(i%) & ", " & STR$(%True) & ", " & STR$(Spec(1,i%)) INCR i% LOOP UNTIL i% = 17 PRINT #nr, $Spectrum_End PRINT #nr, SPACE$(10) PRINT #nr,"<";f;">-file created by on "; DATE$ PRINT #nr,"GMT is written by ";$GWR PRINT #nr,"assisted by Kristof Lauwers - Logos Foundation" PRINT #nr, $EMAIL PRINT #nr, $WEBSITE PRINT #nr,"[EOF]" CLOSE #nr FUNCTION = %True END FUNCTION FUNCTION CreateUniqueFilename (BYVAL appid AS DWORD) EXPORT AS STRING ' 16.01.2001 ' the string returned starts with "g_" so that we can recognize the origin later on. ' This function was written for ftp interactive applications. ' appid should identify the application creating the filename. LOCAL datum AS STRING * 8 LOCAL tijd AS STRING * 6 LOCAL pin AS STRING * 4 LOCAL mpName AS ASCIIZ * %MAX_COMPUTERNAME_LENGTH + 1 GetComputerName mpName,SIZEOF(mpName) ' WinApi function mpname = MCASE$(TRIM$(mpname)) datum = RIGHT$(DATE$,4) & LEFT$(DATE$,2) & MID$(DATE$,4,2) tijd = LEFT$(TIME$,2) & MID$(TIME$,4,2) & RIGHT$(TIME$,2) RANDOMIZE TIMER pin = HEX$(RND(1) * &H0FFFF, 4) FUNCTION = "g_" & LTrimZero(HEX$(appid)) & "_" & datum & "_" & tijd & "_" & pin & "_" & mpname & ".asc" END FUNCTION FUNCTION IsPrime (BYVAL n AS DWORD) EXPORT AS DWORD LOCAL i AS DWORD FUNCTION = %True SELECT CASE n CASE 0,1,2,3,5,7,11,13,17,19,23,29 EXIT FUNCTION END SELECT FOR i = 2 TO n-1 IF ISFALSE n MOD i THEN FUNCTION = %False : EXIT FUNCTION NEXT i END FUNCTION FUNCTION NextPrime (BYVAL n AS DWORD) EXPORT AS DWORD LOCAL retval AS DWORD 'special case for 2 added by xof 30/01/2013 - before nextprime(1) would return 3... IF n = 1 THEN FUNCTION = 2 EXIT FUNCTION END IF IF (n AND 1) THEN retval = n + 2 ELSE retval = n + 1 DO IF IsPrime (retval) THEN EXIT LOOP ELSE retval += 2 END IF LOOP FUNCTION = retval END FUNCTION FUNCTION PrevPrime (BYVAL n AS DWORD) EXPORT AS DWORD LOCAL retval AS DWORD SELECT CASE n CASE 0 FUNCTION = 0 : EXIT FUNCTION ' illegal CASE 1 FUNCTION = 1 : EXIT FUNCTION CASE 2 FUNCTION = 1 : EXIT FUNCTION CASE 3 FUNCTION = 2 : EXIT FUNCTION CASE 5 FUNCTION = 3 : EXIT FUNCTION CASE 7 FUNCTION = 5 : EXIT FUNCTION END SELECT IF (n AND 1) THEN retval = n - 2 ELSE retval = n - 1 DO IF IsPrime (retval) THEN EXIT LOOP ELSE retval -= 2 IF retval = 7 THEN EXIT LOOP END IF END IF LOOP FUNCTION = retval END FUNCTION FUNCTION bitcount(BYVAL x AS DWORD) EXPORT AS LONG ' returns the number of bits that are set in a dword/long LOCAL cnt AS DWORD LOCAL i AS DWORD cnt = %false FOR i = 0 TO 31 IF BIT(x, i) THEN INCR cnt NEXT FUNCTION = cnt END FUNCTION FUNCTION FlipByte (BYVAL b AS BYTE) EXPORT AS BYTE ' convert little endian into big endian numbers LOCAL retval AS BYTE IF BIT (b,0) THEN BIT SET retval,7 IF BIT (b,1) THEN BIT SET retval,6 IF BIT (b,2) THEN BIT SET retval,5 IF BIT (b,3) THEN BIT SET retval,4 IF BIT (b,4) THEN BIT SET retval,3 IF BIT (b,5) THEN BIT SET retval,2 IF BIT (b,6) THEN BIT SET retval,1 IF BIT (b,7) THEN BIT SET retval,0 FUNCTION = retval END FUNCTION FUNCTION FlipWord (BYVAL w AS WORD) EXPORT AS WORD LOCAL retval AS WORD LOCAL b AS BYTE #IF %PB_REVISION < &H800 'then b = FlipByte(LOBYT(w)) #ELSE b = FlipByte(LO(BYTE,w)) #ENDIF retval = b SHIFT LEFT retval,8 #IF %PB_REVISION < &H800 'then b = FlipByte(HIBYT(w)) #ELSE b = FlipByte(HI(BYTE,w)) #ENDIF FUNCTION = retval OR b END FUNCTION FUNCTION FlipDword (BYVAL dwr AS DWORD) EXPORT AS DWORD LOCAL retval AS DWORD LOCAL b AS WORD #IF %PB_REVISION < &H800 'then b = FlipWord(LOWRD(dwr)) #ELSE b = Flipword(LO(WORD,dwr)) #ENDIF retval = b SHIFT LEFT retval,16 #IF %PB_REVISION < &H800 b = FlipWord(HIWRD(dwr)) #ELSE b = Flipword(HI(WORD, dwr)) #ENDIF FUNCTION = retval OR b END FUNCTION FUNCTION FlipQuad (BYVAL q AS QUAD) EXPORT AS QUAD LOCAL retval AS QUAD LOCAL b AS DWORD b = FlipdWord(q AND &HFFFFFFFF&) retval = b SHIFT LEFT retval,32 SHIFT RIGHT q,32 b = FlipDword(q AND &HFFFFFFFF&) FUNCTION = retval OR b END FUNCTION SUB RotateLowNibble (BYREF b AS BYTE, BYVAL direction AS LONG) EXPORT ' 15.09.2002 ' used for stepping motor controls ' rotates the bit pattern in the low nibble of the byte passed a single bit.` ' preserves other bits in the byte passed LOCAL n AS BYTE n = b AND &H0F0 ' save the high nibble in b b AND= &H0F ' isolate the low nibble SELECT CASE direction CASE < 0 ' rotate left ROTATE LEFT b,1 IF BIT (b,4) THEN BIT SET b,0 CASE ELSE ' rotate right ROTATE RIGHT b, 1 IF BIT (b,7) THEN BIT SET b,3 END SELECT b AND= &H0F ' mask high nibble b OR= n ' place the existing high nibble back END SUB SUB RotateHighNibble (BYREF b AS BYTE, BYVAL direction AS LONG) EXPORT ' 15.09.2002 ' used for stepping motor controls ' rotates the bit pattern in the high nibble of the byte passed a single bit.` ' preserves other bits in the byte passed LOCAL n AS BYTE n = b AND &H0F ' save the low nibble in b b = b AND &H0F0 ' isolate the high nibble SELECT CASE direction CASE < 0 ' rotate left ROTATE LEFT b,1 IF BIT (b,0) THEN BIT SET b,4 CASE ELSE ' rotate right ROTATE RIGHT b, 1 IF BIT (b,3) THEN BIT SET b,7 END SELECT b = b AND &H0F0 ' mask low nibble b = b OR n ' place the existing low nibble back END SUB FUNCTION NextE12 (BYVAL i AS SINGLE) EXPORT AS SINGLE LOCAL n AS SINGLE LOCAL decexp AS INTEGER ' returns next higher value in the E12 series ' 1.0, 1.2, 1.5, 1.8, 2.2, 2.7, 3.3, 3.9, 4.7, 5.6, 6.8, 8.2 IF i = %False THEN FUNCTION = %False : EXIT FUNCTION decexp = %False IF i >= 10 THEN DO i /= 10! INCR decexp n = ROUND(i,1) LOOP UNTIL n < 10 ELSEIF i >= 1 THEN decexp = 0 n = ROUND(i,1) ELSE DO i *= 10! DECR decexp n = ROUND(i,1) LOOP UNTIL n >= 1 END IF ' now we should have the decimal exponent in decexp and the mantissa in n SELECT CASE n CASE < 1.2 FUNCTION = 1.2 * (10^decexp) CASE < 1.5 FUNCTION = 1.5 * (10^decexp) CASE < 1.8 FUNCTION = 1.8 * (10^decexp) CASE < 2.2 FUNCTION = 2.2 * (10^decexp) CASE < 2.7 FUNCTION = 2.7 * (10^decexp) CASE < 3.3 FUNCTION = 3.3 * (10^decexp) CASE < 3.9 FUNCTION = 3.9 * (10^decexp) CASE < 4.7 FUNCTION = 4.7 * (10^decexp) CASE < 5.6 FUNCTION = 5.6 * (10^decexp) CASE < 6.8 FUNCTION = 6.8 * (10^decexp) CASE < 8.2 FUNCTION = 8.2 * (10^decexp) CASE < 10 FUNCTION = 10 * (10^decexp) END SELECT END FUNCTION FUNCTION PrevE12 (BYVAL i AS SINGLE) EXPORT AS SINGLE LOCAL n AS SINGLE LOCAL decexp AS INTEGER ' returns next lower value in the E12 series ' 1.0, 1.2, 1.5, 1.8, 2.2, 2.7, 3.3, 3.9, 4.7, 5.6, 6.8, 8.2 IF i = %False THEN FUNCTION = %False : EXIT FUNCTION decexp = %False IF i >= 10 THEN DO i /= 10! INCR decexp n = ROUND(i,1) LOOP UNTIL n < 10 ELSEIF i >= 1 THEN decexp = 0 n = ROUND(i,1) ELSE DO i *= 10! DECR decexp n = ROUND(i,1) LOOP UNTIL n >= 1 END IF ' now we should have the decimal exponent in decexp and the mantissa in n SELECT CASE n CASE > 8.2 FUNCTION = 10 * (10^decexp) CASE > 6.8 FUNCTION = 8.2 * (10^decexp) CASE > 5.6 FUNCTION = 6.8 * (10^decexp) CASE > 4.7 FUNCTION = 5.6 * (10^decexp) CASE > 3.9 FUNCTION = 4.7 * (10^decexp) CASE > 3.3 FUNCTION = 3.9 * (10^decexp) CASE > 2.7 FUNCTION = 3.3 * (10^decexp) CASE > 2.2 FUNCTION = 2.7 * (10^decexp) CASE > 1.8 FUNCTION = 2.2 * (10^decexp) CASE > 1.5 FUNCTION = 1.8 * (10^decexp) CASE > 1.2 FUNCTION = 1.5 * (10^decexp) CASE > 1 FUNCTION = 1.2 * (10^decexp) END SELECT END FUNCTION FUNCTION GetWindir () EXPORT AS STRING ' returns the complete path and name where windows is installed on this machine LOCAL n AS LONG LOCAL i AS DWORD STATIC dummy$ ' method using environ$: IF dummy$ = "" THEN DO INCR i dummy$ = ENVIRON$(i) IF UCASE$(LEFT$(dummy,6)) = "WINDIR" THEN EXIT LOOP END IF LOOP WHILE dummy$ <> "" IF dummy$ = "" THEN Warning "Windir not found in environment string in " + FUNCNAME$, 10000 FUNCTION = "" EXIT FUNCTION ELSE i = LEN(dummy$) dummy$ = RIGHT$(dummy$,i-7) END IF FUNCTION = dummy$ ELSE FUNCTION = dummy$ END IF END FUNCTION FUNCTION GetGmtDir () EXPORT AS STRING '09.06.2002- 11.06.2002 ' returns the complete path and name where gmt is installed on this machine ' updates the default ini file if appropriate. LOCAL n AS LONG LOCAL i AS DWORD LOCAL file AS STRING LOCAL dum$ STATIC dummy$ IF dummy$ = "" THEN ' first we try retrieving this directory in the inifile: file = IniFileName IF Existfile (file) THEN 'msgbox "reading file" n = FREEFILE OPEN file FOR INPUT LOCK SHARED AS #n WHILE NOT EOF(n) INPUT #n, dum$ dum$ = TRIM$(UCASE$(dum$)) IF dum$ = "[GMT_DIR]" THEN INPUT #n, dummy$ FUNCTION = dummy$ CLOSE #n EXIT FUNCTION END IF WEND CLOSE #n dummy$ = "" END IF IF dummy$ = "" THEN ' if not found in the ini file, we search it using ' shelled DOS commands: n = FREEFILE SHELL ENVIRON$("COMSPEC") + " /C DIR/AD > gmtdir.txt" OPEN "gmtdir.txt" FOR INPUT LOCK SHARED AS #n DO LINE INPUT #n, dummy$ dummy$=LTRIM$(dummy$) IF LEFT$(UCASE$(dummy$),9) = "DIRECTORY" THEN EXIT LOOP LOOP UNTIL EOF(n) i = LEN(dummy$) dummy$= RIGHT$(dummy$, i-13) dummy$ = TRIM$(dummy$) FUNCTION = dummy$ CLOSE #n KILL "gmtdir.txt" ' save this in the ini-file: ' write the directory where gmt is installed to the ini file IF existfile(file) THEN n = FREEFILE OPEN file FOR APPEND LOCK SHARED AS #n PRINT #n, " " PRINT #n,"[GMT_DIR],", dummy$ CLOSE #n END IF END IF ELSE FUNCTION = dummy$ END IF END FUNCTION FUNCTION TriangleMath (Tri AS triangletype) EXPORT AS DWORD ' not yet fully debugged ! (01.09.2002) ' precision limited to 0.00001 ' returns msgbox for unacceptable input. LOCAL atemp AS SINGLE FUNCTION = %False ' Tri.a ' hoek tegenover de zijde sa ' Tri.b ' hoek tegenover de zijde sb ' Tri.c ' hoek tegenover de zijde sc ' Tri.sa ' Tri.sb ' Tri.sc ' Tri.ha ' Tri.hb ' Tri.hc ' Tri.S ' surface ' should solve any triangle, given 3 elements. IF (Tri.sa > 0.0) AND (Tri.sb > 0.0) AND (Tri.sc > 0.0) THEN ' 3 zijden gegeven ---> cosinusregel toepassen ' debug o.k. Tri.a = ((Tri.sb * Tri.sb) + (Tri.sc * Tri.sc) - (Tri.sa * Tri.sa)) / (2 * Tri.sb * Tri.sc) Tri.a = Arccos (Tri.a) Tri.b = ((Tri.sa * Tri.sa) + (Tri.sc * Tri.sc) - (Tri.sb * Tri.sb)) / (2 * Tri.sa * Tri.sc) Tri.b = Arccos (Tri.b) Tri.c = ((Tri.sa * Tri.sa) + (Tri.sb * Tri.sb) - (Tri.sc * Tri.sc)) / (2 * Tri.sa * Tri.sb) Tri.c = ArcCos (Tri.c) ' atemp = min(Tri.a,Tri.b,Tri.c) ' if atemp < 0.0 then ' msgbox "Cannot solve triangle",, funcname$ ' exit function ' end if IF ABS((Tri.a + Tri.b + Tri.c) - (ATN(1) * 4.0)) > 0.00001 THEN Warning "Cannot solve triangle in " + FUNCNAME$, 10000 EXIT FUNCTION END IF GOTO fillrest ELSEIF (Tri.sa = 0.0) AND (Tri.sb = 0.0) AND (Tri.sc = 0.0) THEN ' in this case we need 3 angles ' in this case we can only return the sides as proportional values, the largest side normalized to 1. ' debug o.k. IF (Tri.a > 0.0) AND (Tri.b > 0.0) THEN Tri.c = (ATN(1)*4.0) - (Tri.a + Tri.b) ELSEIF (Tri.a > 0.0) AND (Tri.c > 0.0) THEN Tri.b = (ATN(1)*4.0) - (Tri.a + Tri.c) ELSEIF (Tri.b > 0.0) AND (Tri.c > 0.0) THEN Tri.a = (ATN(1)*4.0) - (Tri.b + Tri.c) ELSE EXIT FUNCTION ' unsolvable END IF IF ABS((Tri.a + Tri.b + Tri.c) - (ATN(1) * 4.0)) > 0.00001 THEN Warning "Cannot solve triangle in " + FUNCNAME$, 10000 EXIT FUNCTION END IF atemp = MAX(Tri.a,Tri.b,Tri.c) IF Tri.a = atemp THEN Tri.sa = 1! Tri.sb = Tri.sa * SIN(Tri.b) / SIN(Tri.a) Tri.sc = Tri.sa * SIN(Tri.c) / SIN(Tri.a) GOTO fillrest END IF IF Tri.b = atemp THEN Tri.sb = 1! Tri.sa = Tri.sb * SIN(Tri.a) / SIN(Tri.b) Tri.sc = Tri.sb * SIN(Tri.c) / SIN(Tri.b) GOTO fillrest END IF IF Tri.c = atemp THEN Tri.sc = 1! Tri.sa = Tri.sc * SIN(Tri.a) / SIN(Tri.c) Tri.sb = Tri.sc * SIN(Tri.b) / SIN(Tri.c) GOTO fillrest END IF ELSEIF (Tri.sa > 0.0) AND (Tri.sb = 0.0) AND (Tri.sc = 0.0) THEN ' geval 2 hoeken bekend en 1 zijde IF (Tri.a > 0.0) AND (Tri.b > 0.0) THEN Tri.c = (ATN(1)*4.0) - (Tri.a + Tri.b) ELSEIF (Tri.a > 0.0) AND (Tri.c > 0.0) THEN Tri.b = (ATN(1)*4.0) - (Tri.a + Tri.c) ELSEIF (Tri.b > 0.0) AND (Tri.c > 0.0) THEN Tri.a = (ATN(1) *4.0) - (Tri.b + Tri.c) ELSE EXIT FUNCTION ' unsolvable END IF atemp = MAX(Tri.sa,Tri.sb,Tri.sc) IF Tri.sa = atemp THEN Tri.sb = Tri.sa * SIN(Tri.b) / SIN(Tri.a) Tri.sc = Tri.sa * SIN(Tri.c) / SIN(Tri.a) ' o.k. GOTO fillrest END IF IF Tri.sb = atemp THEN Tri.sa = Tri.sb * SIN(Tri.a) / SIN(Tri.b) Tri.sc = Tri.sb * SIN(Tri.c) / SIN(Tri.b) GOTO fillrest END IF IF Tri.sc = atemp THEN Tri.sa = Tri.sc * SIN(Tri.a) / SIN(Tri.c) Tri.sb = Tri.sc * SIN(Tri.b) / SIN(Tri.c) GOTO fillrest END IF ELSEIF (Tri.sb > 0.0) AND (Tri.sa = 0.0) AND (Tri.sc = 0.0) THEN ' geval 2 hoeken bekend en 1 zijde IF (Tri.a > 0.0) AND (Tri.b > 0.0) THEN Tri.c = (ATN(1)*4.0) - (Tri.a + Tri.b) ELSEIF (Tri.a > 0.0) AND (Tri.c > 0.0) THEN Tri.b = (ATN(1)*4.0) - (Tri.a + Tri.c) ELSEIF (Tri.b > 0.0) AND (Tri.c > 0.0) THEN Tri.a = (ATN(1)*4.0) - (Tri.b + Tri.c) ELSE EXIT FUNCTION ' unsolvable END IF atemp = MAX(Tri.sa,Tri.sb,Tri.sc) IF Tri.sa = atemp THEN Tri.sb = Tri.sa * SIN(Tri.b) / SIN(Tri.a) Tri.sc = Tri.sa * SIN(Tri.c) / SIN(Tri.a) GOTO fillrest END IF IF Tri.sb = atemp THEN Tri.sa = Tri.sb * SIN(Tri.a) / SIN(Tri.b) Tri.sc = Tri.sb * SIN(Tri.c) / SIN(Tri.b) GOTO fillrest END IF IF Tri.sc = atemp THEN Tri.sa = Tri.sc * SIN(Tri.a) / SIN(Tri.c) Tri.sb = Tri.sc * SIN(Tri.b) / SIN(Tri.c) GOTO fillrest END IF ELSEIF (Tri.sc > 0.0) AND (Tri.sa = 0.0) AND (Tri.sb = 0.0) THEN ' geval 2 hoeken bekend en 1 zijde IF (Tri.a > 0.0) AND (Tri.b > 0.0) THEN Tri.c = (ATN(1)*4.0) - (Tri.a + Tri.b) ELSEIF (Tri.a > 0.0) AND (Tri.c > 0.0) THEN Tri.b = (ATN(1)*4.0) - (Tri.a + Tri.c) ELSEIF (Tri.b > 0.0) AND (Tri.c > 0.0) THEN Tri.a = (ATN(1)*4.0) - (Tri.b + Tri.c) ELSE EXIT FUNCTION ' unsolvable END IF atemp = MAX(Tri.sa,Tri.sb,Tri.sc) IF Tri.sa = atemp THEN Tri.sb = Tri.sa * SIN(Tri.b) / SIN(Tri.a) Tri.sc = Tri.sa * SIN(Tri.c) / SIN(Tri.a) GOTO fillrest END IF IF Tri.sb = atemp THEN Tri.sa = Tri.sb * SIN(Tri.a) / SIN(Tri.b) Tri.sc = Tri.sb * SIN(Tri.c) / SIN(Tri.b) GOTO fillrest END IF IF Tri.sc = atemp THEN Tri.sa = Tri.sc * SIN(Tri.a) / SIN(Tri.c) Tri.sb = Tri.sc * SIN(Tri.b) / SIN(Tri.c) GOTO fillrest END IF ' case 1 angle and 2 adjacent sides ELSEIF (Tri.a > 0.0) AND (Tri.sb > 0.0) AND (Tri.sc > 0.0) THEN ' cosinusregel voor berekening van Tri.sa: Tri.sa = SQR((Tri.sb * Tri.sb) + (Tri.sc * Tri.sc) - (2 * Tri.sb * Tri.sc * COS(Tri.a))) Tri.b = Tri.sb * SIN(Tri.a) / Tri.sa Tri.b = ArcSin(Tri.b) Tri.c = Tri.sc * SIN(Tri.a) / Tri.sa Tri.c = ArcSin(Tri.c) GOTO fillrest ELSEIF (Tri.b > 0.0) AND (Tri.sa > 0.0) AND (Tri.sc > 0.0) THEN Tri.sb = SQR((Tri.sa * Tri.sa) + (Tri.sc * Tri.sc) - (2 * Tri.sa * Tri.sc * COS(Tri.b))) Tri.a = Tri.sa * SIN(Tri.b) / Tri.sb Tri.a = ArcSin(Tri.a) Tri.c = Tri.sc * SIN(Tri.b) / Tri.sb Tri.c = ArcSin(Tri.c) GOTO fillrest ELSEIF (Tri.c > 0.0) AND (Tri.sa > 0.0) AND (Tri.sb > 0.0) THEN Tri.sc = SQR((Tri.sa * Tri.sa) + (Tri.sb * Tri.sb) - (2 * Tri.sa * Tri.sb * COS(Tri.c))) Tri.b = Tri.sb * SIN(Tri.c) / Tri.sc Tri.b = ArcSin(Tri.b) Tri.a = Tri.sa * SIN(Tri.c) / Tri.sc Tri.a = ArcSin(Tri.a) GOTO fillrest ' to be done: cases 1 angle and 2 sides given... ELSEIF Tri.a = 0.0 THEN ' in this case we always need Tri.sa IF Tri.sa > 0.0 THEN IF (Tri.b > 0.0) AND (Tri.sb > 0.0) THEN Tri.a = Tri.sa * SIN(Tri.b) / Tri.sb Tri.a = ArcSin(Tri.a) Tri.c = (ATN(1) * 4) - (Tri.a + Tri.b) IF Tri.c < 0.00000# THEN Warning "Cannot solve triangle in " + FUNCNAME$, 10000 EXIT FUNCTION END IF Tri.sc = Tri.sb * SIN(Tri.c) / SIN(Tri.b) GOTO fillrest END IF IF (Tri.c > 0.0) AND (Tri.sc > 0.0) THEN Tri.a = Tri.sa * SIN(Tri.c) / Tri.sc Tri.a = ArcSin(Tri.a) Tri.b = (ATN(1)*4.0) - (Tri.a + Tri.c) IF Tri.b < 0.0# THEN Warning "Cannot solve triangle in " + FUNCNAME$, 10000 EXIT FUNCTION END IF Tri.sb = Tri.sa * SIN(Tri.b) / SIN(Tri.a) GOTO fillrest END IF EXIT FUNCTION ELSE EXIT FUNCTION END IF ELSEIF Tri.b = 0.0 THEN ' in this case we always need Tri.sb IF Tri.sb > 0.0 THEN IF (Tri.a > 0.0) AND (tri.sa > 0.0) THEN Tri.b = Tri.sb * SIN(Tri.a) / Tri.sa Tri.b = Arcsin(Tri.b) Tri.c = (ATN(1)*4.0) - (Tri.a + Tri.b) IF Tri.c < 0.0# THEN Warning "Cannot solve triangle in " + FUNCNAME$, 10000 EXIT FUNCTION END IF Tri.sc = Tri.sb * SIN(Tri.c) / SIN(Tri.b) GOTO fillrest END IF IF (Tri.c > 0.0) AND (Tri.sc > 0.0) THEN Tri.b = Tri.sa * SIN(Tri.c) / Tri.sc Tri.b = ArcSin(Tri.b) Tri.a = (ATN(1)*4.0) - (Tri.b + Tri.c) IF Tri.a < 0.0# THEN Warning "Cannot solve triangle in " + FUNCNAME$, 10000 EXIT FUNCTION END IF Tri.sa = Tri.sc * SIN(Tri.a) / SIN(Tri.c) GOTO fillrest END IF EXIT FUNCTION ELSE EXIT FUNCTION END IF ELSEIF tri.c = 0.0 THEN ' in this case we always need Tri.sc IF Tri.sc > 0.0 THEN IF (Tri.a > 0.0) AND (Tri.sa > 0.0) THEN Tri.c = Tri.sc * SIN(Tri.a) / Tri.sa Tri.c = ArcSin(Tri.c) Tri.b = (ATN(1)*4.0) - (Tri.a + Tri.c) IF Tri.b < 0.0# THEN Warning "Cannot solve triangle in " + FUNCNAME$, 10000 EXIT FUNCTION END IF Tri.sb = Tri.sa * SIN(Tri.b) / SIN(Tri.a) GOTO fillrest END IF IF (Tri.b > 0.0) AND (Tri.sb > 0.0) THEN Tri.c = Tri.sc * SIN(Tri.b) / Tri.sb Tri.c = ArcSin(Tri.c) Tri.a = (ATN(1)*4.0) - (Tri.b + Tri.c) IF Tri.a < 0.0# THEN Warning "Cannot solve triangle in " + FUNCNAME$, 10000 EXIT FUNCTION END IF Tri.sa = Tri.sc * SIN(Tri.a) / SIN(Tri.c) GOTO fillrest END IF EXIT FUNCTION ELSE EXIT FUNCTION END IF ELSE ' no solution possible EXIT FUNCTION END IF fillrest: ' bereken Tri.ah, Tri.bh, Tri.ch, Tri.S ' hoogtelijnen vanuit hoeken a,b,c op de zijden sa, sb, sc Tri.S = Tri.sa * Tri.sb * SIN(Tri.c) / 2! ' surface Tri.ch = SIN(Tri.a) * Tri.sb Tri.ah = SIN(Tri.b) * Tri.sc Tri.bh = SIN(Tri.c) * Tri.sa FUNCTION = %True END FUNCTION FUNCTION ArcSin (BYVAL value AS SINGLE) EXPORT AS SINGLE ' value = sine of an angle ' the function returns the angle in rad. FUNCTION = ATN(Value / SQR(1 - (Value * Value))) END FUNCTION FUNCTION ArcCos (BYVAL value AS SINGLE) EXPORT AS SINGLE FUNCTION = (ATN(1)*4.0) / 2.0 - ATN(Value / SQR(1 - (Value * Value))) END FUNCTION ' muldiv32.h ' Description: ' math routines for 32 BIT SIGNED and unsigned numbers. ' ' MulDiv32(a,b,c) = (a * b) / c (ROUND DOWN, SIGNED) ' MulDivRD(a,b,c) = (a * b) / c (ROUND DOWN, unsigned) ' MulDivRN(a,b,c) = (a * b + c/2) / c (ROUND nearest, unsigned) ' MulDivRU(a,b,c) = (a * b + c-1) / c (ROUND UP, unsigned) ' FUNCTION MulDiv32(BYVAL a AS LONG,BYVAL b AS LONG,BYVAL c AS LONG) EXPORT AS LONG LOCAL answer AS LONG !push eax !push ebx !push ecx !push edx !mov eax, a ; correct. checked. !mov ebx, b !mov ecx, c !imul ebx !idiv ecx !shld edx, eax, 16 ; could'nt we do: shld answer, eax, 16 ??? !mov answer, edx !pop edx !pop ecx !pop ebx !pop eax FUNCTION = answer END FUNCTION FUNCTION MulDivRN(BYVAL a AS DWORD ,BYVAL b AS DWORD,BYVAL c AS DWORD) EXPORT AS DWORD LOCAL answer AS DWORD !push eax !push ebx !push ecx !push edx ' !mov eax,DWORD PTR a ;// IF passed byref. !mov eax, a ' !mov ebx,DWORD PTR b ;// IF b passed BYREF !mov ebx, b ' !mov ecx,DWORD PTR c ;// IF c passed BYREF !mov ecx, c !mul ebx !mov ebx,ecx !shr ebx,1 ;// sar ebx,1 !add eax,ebx !adc edx,0 !div ecx !shld edx, eax, 16 !mov answer, edx !pop edx !pop ecx !pop ebx !pop eax FUNCTION = answer END FUNCTION FUNCTION MulDivRU(BYVAL a AS DWORD,BYVAL b AS DWORD,BYVAL c AS DWORD) EXPORT AS DWORD LOCAL answer AS DWORD !push eax !push ebx !push ecx !push edx 'ASM mov eax,DWORD PTR a ;// !mov eax, a 'ASM mov ebx,DWORD PTR b ;// !mov ebx, b 'ASM mov ecx,DWORD PTR c ;// !mov ecx, c 'ASM mul ebx ;// !mul ebx 'ASM mov ebx,ecx ;// !mov ebx,ecx 'ASM dec ebx ;// !dec ebx 'ASM ADD eax,ebx ;// !add eax,ebx 'ASM adc edx,0 ;// !adc edx,0 'ASM div ecx ;// !div ecx 'ASM shld edx, eax, 16 ;// !shld edx, eax, 16 !mov answer, edx !pop edx !pop ecx !pop ebx !pop eax FUNCTION = answer END FUNCTION FUNCTION MulDivRD(BYVAL a AS DWORD,BYVAL b AS DWORD,BYVAL c AS DWORD) EXPORT AS DWORD LOCAL answer AS DWORD !push eax !push ebx !push ecx !push edx 'ASM mov eax,DWORD PTR a ;// !mov eax, a 'ASM mov ebx,DWORD PTR b ;// !mov ebx, b 'ASM mov ecx,DWORD PTR c ;// !mov ecx, c 'ASM mul ebx ;// !mul ebx 'ASM div ecx ;// !div ecx 'ASM shld edx, eax, 16 ;// !shld edx, eax, 16 !mov answer, edx !pop edx !pop ecx !pop ebx !pop eax FUNCTION = answer END FUNCTION FUNCTION Valtijd (BYVAL hoogte AS SINGLE) EXPORT AS SINGLE ' returns value in seconds, for afstand given in meters. ' valversnelling = 9.81 m/s2 FUNCTION = SQR(hoogte / 9.81) END FUNCTION FUNCTION Valhoogte (BYVAL tijd AS SINGLE) EXPORT AS SINGLE FUNCTION = tijd * tijd * 9.81 END FUNCTION FUNCTION DottedIP$(ip AS LONG) EXPORT ' convert long integer IP number to xxx.xxx.xxx.xxx format LOCAL x AS BYTE PTR x = VARPTR(ip) FUNCTION = FORMAT$(@x) & "." & FORMAT$(@x[1]) & "." & FORMAT$(@x[2]) & "." & FORMAT$(@x[3]) END FUNCTION FUNCTION IpString2nr (ip AS STRING) EXPORT AS LONG ' convert dotted format IP number to long integer format. ' tested 21.11.2002 - o.k. DIM x(3) AS LOCAL BYTE 'dword LOCAL s$ LOCAL n$ LOCAL cnt AS DWORD LOCAL cn AS DWORD LOCAL ret AS DWORD ip = TRIM$(ip) cnt = %False cn = 1 n$ = "" DO s$ = MID$(ip,cn,1) INCR cn SELECT CASE s$ CASE "0" TO "9" n$ = n$ & s$ CASE "." x(cnt) = VAL (n$) n$ = "" INCR cnt CASE ELSE ' error Warning "format error in IP in " + FUNCNAME$, 10000 n$= "" FUNCTION = %False EXIT FUNCTION END SELECT LOOP UNTIL cn > LEN(ip) x(3)= VAL(n$) ret = x(3) SHIFT LEFT ret,8 ret = ret OR x(2) SHIFT LEFT ret, 8 ret = ret OR x(1) SHIFT LEFT ret, 8 ret = ret OR x(0) FUNCTION = ret END FUNCTION FUNCTION RoundFreq (BYVAL freq AS SINGLE) EXPORT AS SINGLE ' since timing resolution is limited to 1 ms under Windows, it makes ' no sense to specify tasks with frequencies such as 780 or 310. ' The first one will be rounded to 1000, the latter to 333 etc... ' This function will return the exact frequency the PC will use. IF freq > 1000.0 THEN freq = 1000.0 ' 1 ms limit IF freq < 0.00027 THEN freq = 0.00027 ' ca. 1 hour FUNCTION = 1000.0 / (1000 \ freq) END FUNCTION FUNCTION Mi2Name (BYVAL n AS DWORD, BYVAL param AS DWORD) EXPORT AS STRING n = n MOD 12 IF ISFALSE param THEN ' italian/french system SELECT CASE n CASE 0 FUNCTION = "Do " CASE 1 FUNCTION = "Do#" CASE 2 FUNCTION = "Re " CASE 3 FUNCTION = "Mib" CASE 4 FUNCTION = "Mi " CASE 5 FUNCTION = "Fa " CASE 6 FUNCTION = "Fa#" CASE 7 FUNCTION = "Sol" CASE 8 FUNCTION = "Lab" CASE 9 FUNCTION = "La " CASE 10 FUNCTION = "Sib" CASE 11 FUNCTION = "Si " END SELECT ELSE SELECT CASE n CASE 0 FUNCTION = "C " CASE 1 FUNCTION = "Cis" CASE 2 FUNCTION = "D " CASE 3 FUNCTION = "Es " CASE 4 FUNCTION = "E " CASE 5 FUNCTION = "F " CASE 6 FUNCTION = "Fis" CASE 7 FUNCTION = "G " CASE 8 FUNCTION = "Gis" CASE 9 FUNCTION = "A " CASE 10 FUNCTION = "Bes" CASE 11 FUNCTION = "B " END SELECT END IF END FUNCTION FUNCTION nrBitChanges (BYVAL a AS DWORD, BYVAL b AS DWORD) EXPORT AS DWORD '19.08.2003 ' counts the number of bits that change in the transition between two numbers LOCAL cnt AS DWORD LOCAL i AS DWORD FOR i = 0 TO 31 IF BIT(a,i) <> BIT(b,i) THEN INCR cnt NEXT i FUNCTION = cnt END FUNCTION FUNCTION Sire_Velo2MidiNoot_old (BYVAL velo AS INTEGER) EXPORT AS SINGLE ' transfer function calculated with Gausfit. gwr.20.04.2005 ' new version in g_midi.inc !!! IF velo < 6 THEN FUNCTION = %False : EXIT FUNCTION velo = MIN(velo,127) ' eerste versie met 5 meetpunten, via Gausfit: 'FUNCTION = 41.12872 + (1.560877 * velo) - (1.880568E-02 * (velo^2!)) + (7.231988E-05 * (velo^3!)) ' 3e graads vergelijking op grond van 15 meetpunten: [29.04.2005] FUNCTION = 44.44099 + (1.31992 * velo) - (1.604728E-02 * (velo^2!)) + (6.405674E-05 * (velo^3!)) END FUNCTION FUNCTION Sire_MidiNoot2Velo_old (BYVAL noot AS SINGLE) EXPORT AS WORD ' new version in g_midi.inc !!! ' will return the 7-bit velo value to send in HIBYT(velo) and the 7-bit LSB in LOBYT(velo) 'if noot > 84 then exit function LOCAL velo AS SINGLE LOCAL retval AS WORD LOCAL f AS SINGLE LOCAL lsb AS WORD IF noot < 48 THEN FUNCTION = %False :EXIT FUNCTION IF noot > 84 THEN FUNCTION = %False: EXIT FUNCTION ' eerste opmeting met 4 meetpunten: ' velo = -1170.342 + (61.77662 * noot) - (1.077227 * (noot ^2)) + (6.265741E-03 * (noot^3)) ' 3e graads vergelijking op grond van 15 meetpunten: [29.04.2005] velo = -1678.095 + (85.66193 * noot) - (1.443499 * (noot ^2)) + (8.094788E-03 * (noot^3)) retval = FIX(velo) ' must have a bug... f = FRAC(velo) * 128 ' now 0-127 SHIFT LEFT retval, 8 lsb = FIX(f) FUNCTION = retval OR lsb END FUNCTION FUNCTION GetNoteNrFromKloktype(BYVAL pVaccaNotes AS kloktype PTR, lb AS DWORD, ub AS DWORD, BYVAL note AS SINGLE, BYVAL tol AS SINGLE) EXPORT AS BYTE 'pass the vaccanotes() array BY POINTER!!! passing it byref caused crashes.. 'actually in this form this function can be used for all instruments that have a kloktype '20060809 changed by xof so that it gives closest match in stead of lowest.. LOCAL nt AS SINGLE LOCAL i AS LONG LOCAL nr AS BYTE LOCAL notes() AS BYTE LOCAL devi() AS SINGLE LOCAL pV() AS Kloktype PTR DIM notes(127) DIM devi(127) DIM pv(127) MAT devi()= CON(100) nt = tol / 100! ' normalize to 0-1 FOR i = lb TO ub notes(i) = i IF (@pVaccaNotes.nf > 0) THEN devi(i) = ABS(@pVaccaNotes.nf - note) pv(i) = pVAccaNotes END IF INCR pVaccaNotes NEXT ARRAY SORT devi(), TAGARRAY notes(), ASCEND IF devi(0) < nt THEN FUNCTION = notes(0) END IF END FUNCTION FUNCTION Scalingfactor_Exp (BYVAL bereik AS SINGLE, BYVAL ex AS SINGLE) EXPORT AS SINGLE ' this function returns the scaling factor such that a function of the form ' x = y ^ (ex) ' returns the same range for x-output values as for the y-input values. ' Reminder: for values of ex > x the curve is increasing difference for higher values of y (cfr. powers) ' for values of ex < 1 the curve is decreasing (cfr roots) ' calculate the function for the maximum value in range: 'x = range ^ ex ' now calculate the rescaling factor: 'function = range / x ' in one step: FUNCTION = bereik / (bereik ^ ex) ' now, calling a the result of this function, the function x = (y ^ ex) * a , will result in the same range for x as for y. END FUNCTION SUB Warning (m AS STRING, OPT BYVAL t AS DWORD) EXPORT ' we create a selfdestructing static window to display the warning. ' m is the text to be displayed ' t is the duration is milliseconds the window should stay on the screen. ' note that WarnTime as well as the message itself are declared as global. 'update 050419: t became optional and defaults to 10000 ' extra dialog doevents to force redrawing 'update 060809: text truncated @ 20 lines STATIC hW AS LONG LOCAL rows AS LONG LOCAL cols AS LONG LOCAL i AS LONG IF ISFALSE t THEN t = 10000 IF ISFALSE WarnTime THEN 'there is currently no warning window visible.. WarnTime = t ' global WarnText = TRIM$(m) + CHR$(0) ' global - used to pass it to the callback. rows = PARSECOUNT(WarnText, $CRLF) FOR i = 1 TO rows: cols = MAX(cols, LEN(PARSE$(TRIM$(WarnText), $CRLF, i))): NEXT cols = (cols + 1) * 4: rows = (rows + 1) * 8 ' if we have already a warning window on screen, we can kill it... IF hW THEN DIALOG END hw DIALOG DOEVENTS END IF DIALOG FONT "Lucida Console", 10 'fixed width DIALOG NEW 0, "warning", 1, 1, cols, rows + 1, %WS_DLGFRAME OR %WS_POPUP OR %DS_MODALFRAME OR %DS_SETFONT OR %WS_BORDER OR %WS_CAPTION OR %WS_SYSMENU TO hw CONTROL ADD TEXTBOX, hw, 1, "" + warntext, 1, 1, cols, rows, %ES_MULTILINE OR %ES_READONLY DIALOG SET COLOR hw, &H0000FF, &H000000 CONTROL SET COLOR hw, 1, &H0000FF, &H000000 'the following line control is here only so we can give it the focus instead of the textbox '(so the text is not selected) ' CONTROL ADD LINE, hw, 2, "", 1, rows + 1, cols, 1, %SS_BLACKFRAME DIALOG SHOW MODELESS hw, CALL WarnProc 'give focus to the line ' CONTROL SET FOCUS hw, 1 ' SendMessage hw, %WM_USER + 1, 0, 0 'creates first timer ELSE 'there is allready an existing window - update it and send msg WM_USER + 1 WarnTime = t WarnText = REMOVE$(WarnText, CHR$(0)) + $CRLF + TRIM$(m) + CHR$(0) rows = PARSECOUNT(WarnText, $CRLF) IF rows > 20 THEN WarnText = MID$(WarnText, INSTR(WarnText, $CRLF) + 2) END IF FOR i = 1 TO rows: cols = MAX(cols, LEN(PARSE$(WarnText, $CRLF, i))): NEXT cols = (cols + 1) * 4: rows = (rows + 2) * 8 DIALOG SET SIZE hw, cols, rows CONTROL SET SIZE hw, 1, cols, rows CONTROL SET TEXT hw, 1, WarnText DIALOG DOEVENTS ' SendMessage hw, %WM_USER + 2, 0, 0 'renews the timer END IF DIALOG DOEVENTS: DIALOG DOEVENTS 'make sure the box is refreshed immediately END SUB CALLBACK FUNCTION WarnProc AS LONG ' not exported !!! STATIC hTimer AS LONG SELECT CASE CBMSG CASE %WM_USER + 1 '%WM_CREATE hTimer = SetTimer(CBHNDL, BYVAL &H0000FEED, BYVAL WarnTime, BYVAL %Null) FUNCTION = %False EXIT FUNCTION CASE %WM_USER + 2 IF hTimer THEN KillTimer %Null, hTimer hTimer = SetTimer(CBHNDL, BYVAL &H0000FEED, BYVAL WarnTime, BYVAL %Null) FUNCTION = %False EXIT FUNCTION CASE %WM_TIMER DIALOG END CBHNDL WarnText = "" WarnTime = 0 CASE %WM_DESTROY IF hTimer THEN KillTimer %Null, hTimer hTimer = 0 END IF WarnText = "" WarnTime = 0 FUNCTION = %False EXIT FUNCTION END SELECT END FUNCTION FUNCTION Velo2Pulse (BYREF robot AS musician, BYVAL velo AS BYTE) EXPORT AS DWORD ' returns value in microseconds based on the ' mapping used by the PIC or GMT for the machine passed. ' This function is required to calculate the maximum possible repetition rate for the ' individual beaters/notes SELECT CASE UCASE$(TRIM$((robot.naam))) CASE "VIBI" FUNCTION = 2048 + (velo * 128 / 4) + ((velo^2) / 4) ' cfr. PIC lookup CASE "TUBI" FUNCTION = 3582 + (velo * 64) + ((velo^2) / 2) END SELECT END FUNCTION FUNCTION LTrimZero (a AS STRING) EXPORT AS STRING ' removes leading zeros from hexadecimal, octal and binary numbers passed as strings ' added 26.03.2005 to overcome changes in the PBWIN compiler. IF LEN(a) < 2 THEN FUNCTION = a : EXIT FUNCTION DO IF LEFT$(a,1) = "0" THEN a = RIGHT$(a,LEN(a)-1) ELSE FUNCTION = a EXIT FUNCTION END IF IF LEN(a) < 2 THEN EXIT LOOP LOOP FUNCTION = a END FUNCTION FUNCTION Bentham (BYVAL n AS SINGLE) EXPORT AS SINGLE ' 12.05.2005 IF n = %False THEN MSGBOX "Illegal call for Bentham distribution",, FUNCNAME$ FUNCTION = %False EXIT FUNCTION END IF FUNCTION = LOG(1.0 + (1.0 / n)) END FUNCTION '----- functions for the design of musical robots and instruments.------------------------------------------ '----- ********************************************************** 'FUNCTION FrustumFlow (BYVAL traject AS SINGLE, BYVAL diam AS SINGLE, BYVAL angle AS SINGLE) EXPORT AS SINGLE ' ' function for the calculation of the optimum orifice for a given conical valve ' ' traject = movement traject of the valve (m) ' ' angle = angle of the cone top in degrees ' ' diameter = largest diameter of the valve cone seat (m)| ' ' the value returned is the diameter of the cilindrical tube or orifice ' ' written and debugged by Godfried-Willem Raes, 30.01.2010 ' ' LOCAL b AS DOUBLE ' LOCAL r1 AS SINGLE ' LOCAL r2 AS SINGLE ' LOCAL r3 AS SINGLE ' LOCAL h AS SINGLE ' b = deg2rad#(angle/2) ' r2 = diam /2 ' r3 = r2 - (SIN(b) * COS(b) * traject) ' h = traject * (SIN(b))^2 ' 'S = Pi# * (r2 + r3) * SQR(((r2-r3)^2) + h ^2) ' ' PRINT "Frustum surface:"; S ; " m^2) surface of thew windflow channel!!! ' ' thus ' r1 = SQR( (r2 + r3) * SQR(((r2-r3)^2) + h^2)) ' FUNCTION = r1 * 2 ' return diameter in meters 'END FUNCTION FUNCTION FrustumArea (BYVAL h AS SINGLE, BYVAL diam1 AS SINGLE, BYVAL diam2 AS SINGLE) EXPORT AS SINGLE ' returns the surface area of the side ' if diam1=diam2 the function returns the side surface of a cylinder. LOCAL r1 AS SINGLE LOCAL r2 AS SINGLE LOCAL area AS SINGLE IF diam1 > diam2 THEN SWAP diam1, diam2 ' diam1 must be smallest r1 = diam1 /2 r2 = diam2 /2 area = Pi * (r1 + r2) * SQR(((r2-r1)^2) + (h^2)) ' Pi# *(r1 + r2)* SQR(((r1-r2)^2) + (h^2)) FUNCTION = area END FUNCTION SUB CalcConeTraject (BYVAL orifice AS SINGLE, BYVAL angle AS SINGLE, Conevalve AS conevalveType) EXPORT LOCAL b, r1, s1, g AS SINGLE Conevalve.possible = %True ConeValve.angle = angle Conevalve.orifice = orifice b = deg2rad#(angle/2) r1 = orifice/ 2 s1 = Pi * (r1^2) ' must be equal to frustum surface Conevalve.area = s1 Conevalve.traject = r1 / (TAN(b) - (SIN(b) * COS(b))) g = SIN(b)* Conevalve.traject Conevalve.diam = Conevalve.orifice + (2 * g * COS(b)) ' minimale waarde Conevalve.depth = (Conevalve.diam/2) * (1/TAN(b)) ' now check conditions: IF Conevalve.diam <= Conevalve.orifice THEN Conevalve.possible = %False IF Conevalve.traject > Conevalve.depth THEN Conevalve.possible = %False END SUB SUB CalcConeOrifice (BYVAL traject AS SINGLE, BYVAL angle AS SINGLE, Conevalve AS conevalveType) EXPORT LOCAL b,k,d3,g,m AS SINGLE Conevalve.possible = %True Conevalve.traject = traject Conevalve.angle = angle b = deg2rad#(angle/2) k = traject * TAN(b) d3 = 2 * k g = SIN(b) * Traject m = g * COS(b) Conevalve.orifice = d3 - (2 * m) Conevalve.area = Pi * ((Conevalve.orifice / 2)^2) ' must be equal to frustum surface Conevalve.diam = Conevalve.orifice + (2 * g * COS(b)) ' minimale waarde Conevalve.depth = (Conevalve.diam/2) * (1/TAN(b)) ' now check conditions: IF Conevalve.diam <= Conevalve.orifice THEN Conevalve.possible = %False IF Conevalve.traject > Conevalve.depth THEN Conevalve.possible = %False END SUB SUB CalcConeAngle (BYVAL traject AS SINGLE,BYVAL orifice AS SINGLE, Conevalve AS ConevalveType, OPTIONAL BYVAL RES AS SINGLE) EXPORT LOCAL r1,s,b,g,k,h,p,d3,r3,angle, fault, e AS SINGLE LOCAL aa, k3, k1, sf, d1 AS SINGLE Conevalve.orifice = orifice Conevalve.traject = traject Conevalve.possible = %True fault = 999999999 IF ISFALSE RES THEN RES = 1 r1 = orifice /2 'orifice d1 = orifice s = Pi * (r1^2) 'this is the surface the aperture frustum should have. Conevalve.area = s angle = 180 - RES ' start from flat valve DO b = deg2rad#(angle/2) k = traject * TAN(b) d3 = 2 * k r3 = k h = traject * ((SIN(b))^2) 'e = ABS(r1 - (traject * (tan(b) - sin(b) * cos(b)))) 'stap voor stap berekening van de oppervlakte van de frustum: aa = deg2rad#(180 - angle) ' a' in de tekening k3 = (Pi * d3 * (r3 / SIN(aa/2)))/2 ' opp. van de gehele kegel waaruit het frustum komt k1 = (Pi * d1 * (r1 / SIN(aa/2)))/2 ' opp. van de af te trekken top sf = k3 - k1 ' opp. van het frustum met schuine zijde g e = ABS(s - sf) 'e = ABS(Conevalve.area - ( (r1+r3)* Pi# * sqr(((r3-r1)^2) + (h^2)) )) IF e < fault THEN fault = e Conevalve.angle = angle Conevalve.diam = d3 Conevalve.depth = r3 * (1/ TAN(b)) END IF angle = angle - RES LOOP UNTIL angle <= RES ' CalcConeTraject (orifice, Conevalve.angle, Conevalve) ' now check conditions: IF Conevalve.diam <= orifice THEN Conevalve.possible = %False IF Conevalve.traject > Conevalve.depth THEN Conevalve.possible = %False END SUB FUNCTION CalcFlatValveDiam (BYVAL orifice AS SINGLE, BYVAL traject AS SINGLE) EXPORT AS SINGLE ' returns the minimum required diameter of the valve pallet FUNCTION = orifice + (2* traject) END FUNCTION FUNCTION CalcFlatValveOrifice (BYVAL diam AS SINGLE, BYVAL traject AS SINGLE) EXPORT AS SINGLE ' returns the orifice for a given valve pallet and traject FUNCTION = diam - (2* traject) ' if negative it is impossible. END FUNCTION FUNCTION CalcFlatValveTraject (BYVAL orifice AS SINGLE, BYVAL diam AS SINGLE) EXPORT AS SINGLE ' Tr = ((orifice/2)^2)/ diam = orifice/4 FUNCTION = MIN(orifice/2, (diam - orifice) / 2) ' if negative it is impossible END FUNCTION FUNCTION CalcFlatValveFlow (BYVAL orifice AS SINGLE, BYVAL traject AS SINGLE) EXPORT AS SINGLE ' we assume here that the diameter of the pallet meets the requirement to be >= 3 *(orifice/2) FUNCTION = MIN(Pi * ((orifice/2) ^2), Pi * orifice * traject) END FUNCTION FUNCTION CalcBallValveArea (BYVAL traject AS SINGLE,BYVAL orifice AS SINGLE, BYVAL diam AS SINGLE)EXPORT AS SINGLE 'gwr 05.02.2010 - debugged 06.02.2010 LOCAL m AS SINGLE LOCAL alfa, beta, gamma, theta AS SINGLE LOCAL r1, rc, x, g, h, r2, d2,hf, area, d1 AS SINGLE r1 = orifice/2 rc = diam/2 'tangamma = Traject / r1 ' tg(c) 'gamma = ATN(traject/r1) ' arctangens ' hoek c 'm = traject /SIN(gamma) ' m= SQR((r1^2) + (Traject^2)) 'sin(alfa/2) = orifice / diam x = orifice/ diam ' sinus van a/2 alfa = ATN( x / SQR(1- (x^2))) ' alfa is a/2 h = rc * COS(alfa) ' m = rc - (Tr + h) g = SQR((r1^2) + ((h + Traject)^2)) - rc ' frustum side r2 = r1 * rc / (rc +g) d2 = 2 * r2 ' now we still have to calculate h, the height of the frustum... beta = 2 * arcsin(r2/rc) ' arcsin is in g_indep.dll theta = HalfPi - (beta/2) ' = 90 - beta/2 theta=e in de tekening 'sin(theta) = hf/ g hf = g * SIN(theta) 'area = FrustumArea (hf, d1, d2) ' func. has bug! area = Pi *(r1 + r2)* SQR(((r1-r2)^2) + (hf^2)) FUNCTION = area END FUNCTION FUNCTION Sigmoid (BYVAL lowlim AS SINGLE, BYVAL highlim AS SINGLE, BYVAL value AS SINGLE) EXPORT AS SINGLE ' sigmoid function with autonormalizing, so it always returns values between 0 and 1 LOCAL bereik AS SINGLE LOCAL param AS SINGLE IF lowlim > highlim THEN SWAP lowlim, highlim IF value > highlim THEN FUNCTION = 1 : EXIT FUNCTION IF value < lowlim THEN FUNCTION = 0: EXIT FUNCTION bereik = highlim - lowlim IF lowlim < 0 THEN ' bipolar input data expected, centered around zero param = (value * 12 / bereik) ' normalize input to function to -6 to + 6 ELSE ' unipolar data expected (0 --- positive maximum value) param = 12 * ((value - (bereik/2)) / bereik) ' shift and normalize to -6 ... +6 END IF ' and now the sigmoid transfer function: FUNCTION = 1.0 / (1.0 + ((EXP(1)) ^ (-param))) ' note: e = exp(1) END FUNCTION FUNCTION LogFile(text$, OPT BYVAL filenam$) EXPORT AS LONG 'if called for the first time and/or filenam$ is passed, we open a new file 'date- and timestamp added, extension is allways .log STATIC flog AS LONG STATIC fileopen AS LONG 'init flag - we can't use flog as 0 is a valid filenr IF ISFALSE fileopen THEN fileopen = 1 IF ISFALSE LEN(filenam$) THEN filenam$ = "gmt_debug" END IF IF LEN(filenam$) THEN IF fileopen THEN CLOSE flog filenam$ = PARSE$(filenam$, ".", 1) 'strip extensions REPLACE " " WITH "_" IN filenam$ filenam$ = filenam$ + PARSE$(DATE$, "-", 3) + PARSE$(DATE$, "-", 1) + PARSE$(DATE$, "-", 2) + REMOVE$(TIME$, ANY ":") + ".log" 'adds yyyymmddhhmmss.log flog = FREEFILE OPEN filenam$ FOR OUTPUT ACCESS WRITE LOCK WRITE AS flog END IF PRINT# flog, text$ FUNCTION = flog 'so we can write to it from other places also - good for e.g. dumping an array END FUNCTION FUNCTION BMI (lengte AS SINGLE, gewicht AS SINGLE) EXPORT AS SINGLE ' lengte in meter ' gewicht in kg FUNCTION = gewicht / (lengte^2) END FUNCTION FUNCTION CheckTimerResolution () EXPORT AS DWORD ' we check whether this PC allows us to run the timer with ms resolution... LOCAL tc AS TIMECAPS LOCAL wTimerRes AS WORD LOCAL retval AS LONG LOCAL m AS ASCIIZ * 255 'LOCAL szTitelbox AS ASCIIZ * 30 STATIC tog AS BYTE 'szTitelBox = " timer test..." m = "Hardware timer resolution on this PC" + CHR$(13) retval = timeGetDevCaps (tc,SIZEOF(tc)) ' winApi call wTimerRes = MIN(MAX(tc.wPeriodMin, %TARGET_RESOLUTION), tc.wPeriodMax) IF ISFALSE tog THEN IF retval = %TIMERR_NOCANDO THEN ' win constant m = m + "is insufficient to run jitterfree." + CHR$(13) m = m + "The Win API reports " + STR$(tc.wPeriodMin) + "ms as minimum value" + CHR$(13) 'MessageBox hInst,m, szTitelbox,%MB_OK OR %MB_ICONSTOP OR %MB_TASKMODAL OR %MB_TOPMOST MSGBOX m,,FUNCNAME$ ELSE #IF %DEF(%Wordy) IF %Wordy > 2 THEN m = m + "is sufficient to run properly." + CHR$(13) m = m + "The Win API reports " + STR$(tc.wPeriodMin) + "ms as minimum value" + CHR$(13) 'MessageBox hInst,m, szTitelbox,%MB_OK OR %MB_ICONASTERISK OR %MB_TASKMODAL OR %MB_TOPMOST MSGBOX m,, FUNCNAME$ END IF #ENDIF END IF END IF IF ISFALSE tog THEN tog = %True ' timeBeginPeriod wTimerRes - we leave this to the application. ' If used, must be matched with a call to ' timeEndPeriod wTimerRes FUNCTION = tc.wPeriodMin END FUNCTION FUNCTION WaveFreq (Samp() AS SINGLE, BYVAL samplingrate AS LONG, BYVAL noise AS SINGLE) EXPORT AS SINGLE ' 14.06.2009 - gwr. tested and debugged ' note that the data must be bipolar!!! ' samplingrate: in S/s ' noise: for a 16 bit converter, this should be not smaller than 1/32768 = 3E-5 ' for a 12 bit converter, 1/ 4096 = 2.44 E-4 ' for 10 bit resolution, 1/ 1024 = ca.0.001 ' note that the data array does not have to be normalized to -1 - +1, but ' if it is not, the noise values have to be recalculated. LOCAL time AS SINGLE LOCAL oldsign AS LONG LOCAL i, j, cnt AS DWORD time = (UBOUND(Samp) + 1)/ samplingrate ' tijdsduur van het aangeboden array, in sekonden. DO IF ABS(Samp(i)) > noise THEN oldsign = SGN(Samp(i)) IF oldsign <> %False THEN EXIT DO END IF INCR i LOOP UNTIL i >= UBOUND(samp) IF i >= UBOUND(samp) THEN FUNCTION = %False: EXIT FUNCTION ' only noise in buffer INCR i FOR j = i TO UBOUND(Samp) IF ABS(samp(j)) > noise THEN IF SGN(samp(j)) <> oldsign THEN INCR cnt ' cnt is incremented on each zerocross above trigger oldsign = SGN(samp(j)) END IF ELSE ITERATE FOR END IF NEXT j ' omrekening naar absolute frekwentie: IF cnt THEN INCR cnt ELSE FUNCTION = %False: EXIT FUNCTION 'cnt = cnt / 2 'FUNCTION = cnt / time ' in Hz ' coded a bit faster: FUNCTION = cnt / (time + time) END FUNCTION FUNCTION WaveFreq_Dbl (Samp() AS DOUBLE, BYVAL samplingrate AS LONG, BYVAL noise AS DOUBLE) EXPORT AS SINGLE ' 14.06.2009 - gwr. tested and debugged ' note that the data must be bipolar!!! ' samplingrate: in S/s ' noise: for a 16 bit converter, this should be not smaller than 1/32768 = 3E-5 ' for a 12 bit converter, 1/ 4096 = 2.44 E-4 ' for 10 bit resolution, 1/ 1024 = ca.0.001 ' note that the data array does not have to be normalized to -1 - +1, but ' if it is not, the noise values have to be recalculated. ' This version works with an array of doubles. LOCAL time AS SINGLE LOCAL oldsign AS LONG LOCAL i, j AS DWORD LOCAL cnt AS DWORD time = (UBOUND(Samp) + 1)/ samplingrate ' tijdsduur van het aangeboden array, in sekonden. i = 0 cnt = 0 oldsign = 0 DO IF ABS(Samp(i)) > noise THEN oldsign = SGN(Samp(i)) IF oldsign <> %False THEN EXIT DO END IF INCR i LOOP UNTIL i >= UBOUND(samp) IF i >= UBOUND(samp) THEN FUNCTION = %False: EXIT FUNCTION ' only noise in buffer INCR i FOR j = i TO UBOUND(Samp) IF ABS(samp(j)) > noise THEN IF SGN(samp(j)) <> oldsign THEN INCR cnt ' cnt is incremented on each zerocross above trigger oldsign = SGN(samp(j)) END IF ELSE ITERATE FOR END IF NEXT j ' omrekening naar absolute frekwentie: IF cnt THEN INCR cnt ELSE FUNCTION = %False: EXIT FUNCTION cnt /= 2 FUNCTION = cnt / time ' in Hz END FUNCTION FUNCTION WaveFreq_Int (Samp() AS INTEGER, BYVAL samplingrate AS LONG, BYVAL noise AS WORD) EXPORT AS SINGLE ' 14.06.2009 - gwr. tested and debugged ' note that the data must be bipolar!!! ' samplingrate: in S/s ' noise: for a 16 bit converter, this should be not smaller than 1 ' This version works with an array of integers. LOCAL time AS SINGLE LOCAL oldsign AS LONG LOCAL i, j AS DWORD LOCAL cnt AS DWORD time = (UBOUND(Samp) + 1)/ samplingrate ' tijdsduur van het aangeboden array, in sekonden. i = 0 cnt = 0 oldsign = 0 DO IF ABS(Samp(i)) > noise THEN oldsign = SGN(Samp(i)) IF oldsign <> %False THEN EXIT DO END IF INCR i LOOP UNTIL i >= UBOUND(samp) IF i >= UBOUND(samp) THEN FUNCTION = %False: EXIT FUNCTION ' only noise in buffer INCR i FOR j = i TO UBOUND(Samp) IF ABS(samp(j)) > noise THEN IF SGN(samp(j)) <> oldsign THEN INCR cnt ' cnt is incremented on each zerocross above trigger oldsign = SGN(samp(j)) END IF ELSE ITERATE FOR END IF NEXT j ' omrekening naar absolute frekwentie: IF cnt THEN INCR cnt ELSE FUNCTION = %False: EXIT FUNCTION cnt = cnt / 2 FUNCTION = cnt / time ' in Hz END FUNCTION FUNCTION WaveFreq_Lng (Samp() AS LONG, BYVAL samplingrate AS LONG, BYVAL noise AS DWORD) EXPORT AS SINGLE ' 14.06.2009 - gwr. tested and debugged ' note that the data must be bipolar!!! ' samplingrate: in S/s ' noise: in units, depending on range in Samp() ' This version works with an array of 32-bit integers. LOCAL time AS SINGLE LOCAL oldsign AS LONG LOCAL i, j AS DWORD LOCAL cnt AS DWORD time = (UBOUND(Samp) + 1)/ samplingrate ' tijdsduur van het aangeboden array, in sekonden. i = 0 cnt = 0 oldsign = 0 DO IF ABS(Samp(i)) > noise THEN oldsign = SGN(Samp(i)) IF oldsign <> %False THEN EXIT DO END IF INCR i LOOP UNTIL i >= UBOUND(samp) IF i >= UBOUND(samp) THEN FUNCTION = %False: EXIT FUNCTION ' only noise in buffer INCR i FOR j = i TO UBOUND(Samp) IF ABS(samp(j)) > noise THEN IF SGN(samp(j)) <> oldsign THEN INCR cnt ' cnt is incremented on each zerocross above trigger oldsign = SGN(samp(j)) END IF ELSE ITERATE FOR END IF NEXT j ' omrekening naar absolute frekwentie: IF cnt THEN INCR cnt ELSE FUNCTION = %False: EXIT FUNCTION cnt = cnt / 2 FUNCTION = cnt / time ' in Hz END FUNCTION FUNCTION SortDftDbl (dftarr() AS DOUBLE, binsize AS SINGLE,sp AS Spectrum128DoubleType) EXPORT AS SINGLE ' binsize: voor een 256 punten FFT krijgen we 128 frekwentiebandjes. ' wanneer de sampling rate 256 S/s is, is de binsize = 1Hz , de f-bandjes lopen van 1Hz tot 128Hz ' bij 64 S/s wordt de binsize dan 0.25Hz en de f-bandjes lopen van 0.25Hz tot 32Hz LOCAL i AS DWORD DIM fbins(127) AS SINGLE AT VARPTR (sp.freq(0)) DIM pows(127) AS DOUBLE AT VARPTR (sp.pow(0)) FOR i = 0 TO 127 fbins(i) = (i+1) * binsize ' tagarray pows(i) = dftarr(i) ' copy NEXT i ARRAY SORT pows(),TAGARRAY fbins(), DESCEND ' this sorts sp.freq() and sp.pow() IF sp.pow(0) > (1/128) THEN FUNCTION = sp.freq(0) ' return strongest ELSE FUNCTION = %False END IF END FUNCTION FUNCTION SortDft (dftarr() AS SINGLE, binsize AS SINGLE,sp AS Spectrum128SingleType) EXPORT AS SINGLE ' binsize: voor een 256 punten FFT krijgen we 128 frekwentiebandjes. ' wanneer de sampling rate 256 S/s is, is de binsize = 1Hz , de f-bandjes lopen van 1Hz tot 128Hz ' bij 64 S/s wordt de binsize dan 0.25Hz en de f-bandjes lopen van 0.25Hz tot 32Hz LOCAL i AS DWORD DIM fbins(127) AS SINGLE AT VARPTR (sp.freq(0)) DIM pows(127) AS SINGLE AT VARPTR (sp.pow(0)) FOR i = 0 TO 127 fbins(i) = (i+1) * binsize ' tagarray pows(i) = dftarr(i) ' copy NEXT i ARRAY SORT pows(),TAGARRAY fbins(), DESCEND ' this sorts sp.freq() and sp.pow() IF sp.pow(0) > (1/128) THEN FUNCTION = sp.freq(0) ' return strongest ELSE FUNCTION = %False END IF END FUNCTION FUNCTION OctaveBands_Dbl (Sp() AS DOUBLE, Okt() AS DOUBLE) EXPORT AS DOUBLE '13.05.2010: reduces a spectrum to the power-bands per octave. ' there is no error checking here on the size of Sp() ' used and tested in the gesture.airborne detection code. LOCAL i,j,k, s AS DWORD LOCAL p AS DOUBLE s = LOG2(UBOUND(Sp)+1) 'thus for Sp(127), s will be 7 IF UBOUND(Okt) < s-1 THEN REDIM Okt(s-1) RESET Okt() FOR i = 0 TO s-1 'for a 128 band spectrum, i runs 0 - 6 j = EXP2(i) - 1 'j runs 0,1,3,7,15,31,63 FOR k = j TO EXP2(i+1) - 2 okt(i) += Sp(k) NEXT k p += okt(i) NEXT i FUNCTION = p 'return the sum off all spectral bands END FUNCTION ' array minimum and maximum functions added 13.07.2009 - gwr FUNCTION ArrayMax_SNG (ar() AS SINGLE) EXPORT AS LONG ' this function returns the i pointer for the first highest value in the array LOCAL i, p AS LONG LOCAL mx AS SINGLE mx = -8.43E-37 FOR i = LBOUND(ar) TO UBOUND(ar) IF ar(i) > mx THEN mx = ar(i) : p = i NEXT i FUNCTION = p END FUNCTION FUNCTION ArrayMax_DBL (ar() AS DOUBLE) EXPORT AS LONG ' this function returns the i pointer for the first highest value in the array LOCAL i, p AS LONG LOCAL mx AS DOUBLE mx = -4.19E-307 FOR i = LBOUND(ar) TO UBOUND(ar) IF ar(i) > mx THEN mx = ar(i) : p = i NEXT i FUNCTION = p END FUNCTION FUNCTION ArrayMax_LNG (ar() AS LONG) EXPORT AS LONG ' this function returns the i pointer for the first highest value in the array LOCAL i AS LONG LOCAL mx AS LONG LOCAL p AS LONG mx = -2147483648 FOR i = LBOUND(ar) TO UBOUND(ar) IF ar(i) > mx THEN mx = ar(i) : p = i NEXT i FUNCTION = p END FUNCTION FUNCTION ArrayMax_INT (ar() AS INTEGER) EXPORT AS LONG ' this function returns the i pointer for the first highest value in the array LOCAL i AS LONG LOCAL mx AS INTEGER LOCAL p AS LONG mx = -32768 FOR i = LBOUND(ar) TO UBOUND(ar) IF ar(i) > mx THEN mx = ar(i) : p = i NEXT i FUNCTION = p END FUNCTION FUNCTION ArrayMin_SNG (ar() AS SINGLE) EXPORT AS LONG ' this function returns the i pointer for the first lowest value in the array LOCAL i AS LONG LOCAL mx AS SINGLE LOCAL p AS LONG mx = 3.4E38 FOR i = LBOUND(ar) TO UBOUND(ar) IF ar(i) < mx THEN mx = ar(i) : p = i NEXT i FUNCTION = p END FUNCTION FUNCTION ArrayMin_DBL (ar() AS DOUBLE) EXPORT AS LONG ' this function returns the i pointer for the first lowest value in the array LOCAL i AS LONG LOCAL mx AS DOUBLE LOCAL p AS LONG mx = 1.79E308 FOR i = LBOUND(ar) TO UBOUND(ar) IF ar(i) < mx THEN mx = ar(i) : p = i NEXT i FUNCTION = p END FUNCTION FUNCTION ArrayMin_LNG (ar() AS LONG) EXPORT AS LONG ' this function returns the i pointer for the first lowest value in the array LOCAL i AS LONG LOCAL mx AS LONG LOCAL p AS LONG mx = 2147483647 FOR i = LBOUND(ar) TO UBOUND(ar) IF ar(i) < mx THEN mx = ar(i) : p = i NEXT i FUNCTION = p END FUNCTION FUNCTION ArrayMin_INT (ar() AS INTEGER) EXPORT AS LONG ' this function returns the i pointer for the first lowest value in the array LOCAL i AS LONG LOCAL mx AS INTEGER LOCAL p AS LONG mx = 32767 FOR i = LBOUND(ar) TO UBOUND(ar) IF ar(i) < mx THEN mx = ar(i) : p = i NEXT i FUNCTION = p END FUNCTION ' ---- in the works FUNCTION ArrayReduce_DBL (ar() AS DOUBLE , BYVAL n AS DWORD, aro() AS DOUBLE) EXPORT AS LONG ' returns the array passed with only the n highest values it contains on entry ' to do: the function returns the number of different values returned REDIM aro(LBOUND(ar) TO UBOUND(ar)) AS DOUBLE LOCAL i AS DWORD LOCAL cnt AS DWORD LOCAL p AS LONG DO p = ArrayMax_DBL (ar()) aro(p) = ar(p) ar(p)= -4.19E-307 ' we remove this one from the array passed on input INCR i LOOP UNTIL i = n -1 ' and now, we could place the values back in ar() as an option... END FUNCTION FUNCTION ArrayReduce_SNG (ar() AS SINGLE , BYVAL n AS DWORD, aro() AS SINGLE) EXPORT AS LONG ' returns the array passed with only the n highest values it contains on entry ' to do: the function returns the number of different values returned REDIM aro(LBOUND(ar) TO UBOUND(ar)) AS SINGLE LOCAL i AS DWORD LOCAL cnt AS DWORD LOCAL p AS LONG DO p = ArrayMax_SNG (ar()) aro(p) = ar(p) ar(p)= -8.43E-37 ' we remove this one from the array passed INCR i LOOP UNTIL i = n -1 END FUNCTION FUNCTION RangeLimit (BYVAL minval AS LONG, BYVAL maxval AS LONG, normval AS SINGLE) EXPORT AS LONG ' added by gwr. 14.10.2011 ' normval is expected to be a normalised parameter in the range 0 to 1 LOCAL bereik AS LONG normval = MIN(MAX(normval,0),1) IF minval > maxval THEN SWAP minval, maxval bereik = maxval - minval FUNCTION = minval + (bereik * normval) END FUNCTION FUNCTION Get1Fcoefs (nrvals AS WORD) EXPORT AS DWORD ' 30.03.2012 - gwr. This gives precize integer values. ' returns a pointer to a TwelveWord_Type ' coefficienten voor perfekte 1/f distributies over nrvals elementen ' Read as: in a set of x values, element 1 must occur Tw.a times, element 2 Tw.b times, element 3 Tw.c times etc... STATIC Tw AS TwelveWord_Type 'local i, j as word RESET Tw SELECT CASE nrvals CASE 1 : Tw.a = 1 ' kan eigenlijk niet CASE 2 : Tw.a =2 : Tw.b =1 CASE 3 : Tw.a =6 : Tw.b =3 : Tw.c = 2 CASE 4 : Tw.a =12 : Tw.b =6 : Tw.c = 4 : Tw.d = 3 CASE 5 : Tw.a = 60 : Tw.b = 30 : Tw.c = 20 : Tw.d = 15 : Tw.e = 12 CASE 6 Tw.a = 60 : Tw.b = 30 : Tw.c = 20 : Tw.d = 15 : Tw.e = 12 : Tw.f = 10 ' algoritmisch: 'DIM a(5) as word at varptr(Tw) 'for i = 0 to 5 : a(i) = i+1 : NEXT i 'j = kgv_arw (a()) 'for i = 0 to 5 : a(i)= j / i+1 : next i CASE 7 : Tw.a = 420 : Tw.b = 210 : Tw.c = 140 : Tw.d = 105 : Tw.e = 84 : Tw.f = 70 : Tw.g = 60 CASE 8 : Tw.a = 840 : Tw.b = 420 : Tw.c = 240 : Tw.d = 210 : Tw.e = 168 : Tw.f = 140 : Tw.g = 120 : Tw.h = 105 CASE 9 '9 * 8 * 7 * 5 Tw.a = 2520 : Tw.b = 1260 : Tw.c = 840 : Tw.d = 630 : Tw.e = 504 : Tw.f = 420 : Tw.g = 360 : Tw.h = 315 : Tw.i = 280 CASE 10 Tw.a = 10 * 9 * 4 * 7 ' 2520 '5040 ' we could use the kgv function here. Tw.b = Tw.a/2 ' 1260 '2520 Tw.c = Tw.a/3 ' 840 '1680 Tw.d = Tw.a/4 ' 630 '1260 Tw.e = Tw.a/5 ' 504 '1008 Tw.f = Tw.a/6 ' 420 '840 Tw.g = Tw.a/7 ' 360 '720 Tw.h = Tw.a/8 ' 315 '630 Tw.i = Tw.a/9 ' 280 '560 Tw.j = Tw.a/10 ' 252 '504 CASE 11 Tw.a = 11 * 10 * 9 * 7 * 4 ' = 27720 - fits word Tw.b = Tw.a/2 '13860 Tw.c = Tw.a/3 '9240 Tw.d = Tw.a/4 '6930 Tw.e = Tw.a/5 '5544 Tw.f = Tw.a/6 '4620 Tw.g = Tw.a/7 ' 3960 Tw.h = Tw.a/8 ' 3465 Tw.i = Tw.a/9 ' 3080 Tw.j = Tw.a/10 ' 2772 Tw.k = Tw.a/11 ' 2520 CASE 12 Tw.a = 11 * 10 * 9 * 7 * 4 ' = 27720 - fits word Tw.b = Tw.a/2 '13860 Tw.c = Tw.a/3 '9240 Tw.d = Tw.a/4 '6930 Tw.e = Tw.a/5 '5544 Tw.f = Tw.a/6 '4620 Tw.g = Tw.a/7 ' 3960 Tw.h = Tw.a/8 ' 3465 Tw.i = Tw.a/9 ' 3080 Tw.j = Tw.a/10 ' 2772 Tw.k = Tw.a/11 ' 2520 Tw.l = Tw.a/12 ' 2310 END SELECT FUNCTION = VARPTR(Tw) END FUNCTION FUNCTION InitBarUnits (OPT BYVAL param AS DWORD) EXPORT AS LONG ' Call this one prior to any FindBarUnit function. ' It improves a lot on speed. ' param: for making optional logs ' evt. for limiting the number of initialisations IF param THEN LogFlag = %True ' LogFlag is global here. FindBarUnits2 (0) FindBarUnits3 (0) FindBarUnits4 (0) FindBarUnits5 (0) FindBarUnits6 (0) RESET LogFlag FUNCTION = %True END FUNCTION FUNCTION FindBarUnits2 (BYVAL nrtiks AS LONG) EXPORT AS DWORD ' returns a pointer to a TwelveWord_Type ' find positive integer solutions for the equation '(2a + b) MOD grouping = %False for a given value of nrtiks. ' this maps the 1/f distribution on a set of only two elements. ' The histogram is perfect with these coefficients. ' note that all variables must have a different value for valid solutions. ' Can we create the permutations and bar groupings algorithmically? ' In this case the minimum value for nrtiks seems to be 8, when grouping is 8. ' The simplest example for this metre is the dactylus: short-short-long or long-short-short ' logfile "FindBarUnits2 - 2 metric values " '- grouping in " & STR$(grouping) ' 29.03.2012: debug o.k. ' to do: sort ab pairs in inverse order of likelyness ' 30.03.2012: dif parameter introduced. LOCAL a,b, sum, i , a0, b0 AS LONG STATIC tog, cnt AS INTEGER STATIC dif AS SINGLE STATIC RiHiSize, MaxPermuts, MinTiks AS DWORD ' for debug and monitoring LOCAL pTw AS TwelveWord_Type PTR IF ISFALSE tog THEN ' in this case we calculate the entire lookup DIM RiHi(48) AS STATIC RiHi_Type ' size 48 is enough here. pTw = Get1FCoefs (2) ' in g_indep.dll RiHi(0).coefa = @pTw.a '2 RiHi(0).coefb = @pTw.b '1 dif = 0.25 RESET cnt, RiHiSize FOR a = 1 TO 16 '32 van een zestiende tot een hele noot, dit zijn het aantal tiks waaruit de waarde van a bestaat. FOR b = 1 TO 16 '32 'IF a <> b THEN IF (ABS(a-b) > MAX(a,b) * dif) THEN ' more conditions required: ' a,b if different from 1, should be divisible by 2 or 3 IF a = 1 THEN a0 = %True ELSE RESET a0 IF (ISFALSE (a MOD 2)) OR (ISFALSE (a MOD 3)) THEN a0 = %True END IF IF b = 1 THEN b0 = %True ELSE RESET b0 IF (ISFALSE (b MOD 2)) OR (ISFALSE (b MOD 3)) THEN b0 = %True END IF IF a0 + b0 = 2 THEN 'sum = (2 * a) + b , of in de veralgemeende formulering: sum = (RiHi(0).coefa * a ) + (Rihi(0).coefb * b) ' we kunnen een voorwaarde invoeren om de aantallen te beperken en de indeling in groepen ' te vereenvoudigen: IF (ISFALSE sum MOD 2) OR (ISFALSE sum MOD 3) THEN IF ISFALSE sum MOD FIX(LOG2(sum)-2) THEN ' 02.04.2012 INCR RiHi(sum).aantal ' bevat het aantal maal dat deze ticksum voorkomt IF RiHi(sum).aantal > MaxPermuts THEN MaxPermuts = RiHi(sum).aantal IF sum > RiHiSize THEN RiHiSize = sum RiHi(sum).Tiks(RiHi(sum).aantal-1).a = a RiHi(sum).Tiks(RiHi(sum).aantal-1).b = b END IF END IF END IF END IF NEXT b NEXT a ' clean up to save memory: REDIM PRESERVE RiHi(RiHiSize) AS STATIC Rihi_Type ' for debug: IF LogFlag THEN logfile "possible ticksums and number of occurences for 2-value sets" logfile "Minimum size for RiHi: " & STR$(RiHiSize) logfile "Minimum size for nr. permutations : " & STR$(Maxpermuts) END IF FOR i = 0 TO UBOUND(RiHi) RiHi(i).coefa = RiHi(0).coefa RiHi(i).coefb = RiHi(0).coefb IF RiHi(i).aantal THEN IF ISFALSE MinTiks THEN MinTiks = i IF logflag THEN logfile STR$(i) & " occurs " & STR$(RiHi(i).aantal) & " times" FOR cnt = 0 TO RiHi(i).aantal -1 logfile " " & STR$(RiHi(i).tiks(cnt).a) & STR$(RiHi(i).tiks(cnt).b) NEXT cnt END IF END IF NEXT i tog = %True IF logflag THEN logfile "Function 2 initialized " & FUNCNAME$ logfile "Minimum value for number of ticks : " & STR$(MinTiks) END IF END IF ' normal runtime mode after initialisation: IF nrtiks < = RiHiSize THEN 'UBOUND(RiHi) THEN SELECT CASE RiHi(nrtiks).aantal CASE %False ' check whether the requested number of ticks has a solution for 1/f FUNCTION = VARPTR(RiHi(nrtiks).Tiks(0)) EXIT FUNCTION CASE 1 ' here we should return the only one answer... FUNCTION = VARPTR(RiHi(nrtiks).Tiks(0)) CASE ELSE ' here filters can be applied as well. DO cnt = (RND(1) ^ SQR(2)) * (RiHi(nrtiks).aantal-1) ' 1/f random distribution LOOP UNTIL RiHi(nrtiks).Tiks(cnt).a > 0 FUNCTION = VARPTR(RiHi(nrtiks).Tiks(cnt)) END SELECT ELSE logfile "nrtiks too large for " & FUNCNAME$ warning "Illegal value for nrtiks passed to " & FUNCNAME$ END IF END FUNCTION FUNCTION FindBarUnits3 (BYVAL nrtiks AS LONG) EXPORT AS DWORD ' returns a pointer to a TwelveWords_Type ' find positive integer solutions for the equation '(6a + 3b + 2c) MOD grouping = %False for the minimum value of nrtiks. ' this maps the 1/f distribution on a set of only three elements. ' The histogram is perfect with these coefficients. ' note that all variables must have a different value for valid solutions. ' logfile "FindBarUnits3 - 3 metric values " ' 29.03.2012: debug o.k. ' to do: sort abc groups in inverse order of likelyness LOCAL a,b,c, sum, i , a0, b0,c0 AS LONG STATIC tog, cnt AS INTEGER STATIC dif AS SINGLE STATIC RiHiSize, MaxPerMuts, MinTiks AS DWORD LOCAL pTw AS TwelveWord_Type PTR IF ISFALSE tog THEN ' in this case we calculate the entire lookup DIM RiHi(170) AS STATIC RiHi_Type ' size 169 is minimum required ' RiHi(0).coefa = 6 ' RiHi(0).coefb = 3 ' RiHi(0).coefc = 2 pTw = Get1FCoefs (3) ' in g_indep.dll RiHi(0).coefa = @pTw.a ' = 6 RiHi(0).coefb = @pTw.b '3 RiHi(0).coefc = @pTw.c '2 dif = 0.25 RESET cnt, RiHiSize FOR a = 1 TO 16 '32 van een zestiende tot een hele noot, dit zijn het aantal tiks waaruit de waarde van a bestaat. FOR b = 1 TO 16 '32 FOR c = 1 TO 16 'IF a <> b AND b <> c AND a <> c THEN IF (ABS(a-b) > MAX(a,b) * dif) AND (ABS(a-c) > MAX(a,c) * dif) AND (ABS(b-c) > MAX(b,c) * dif) THEN ' more conditions required: ' a,b if different from 1, should be divisible by 2 or 3 IF a = 1 THEN a0 = %True ELSE RESET a0 IF (ISFALSE (a MOD 2)) OR (ISFALSE (a MOD 3)) THEN a0 = %True END IF IF b = 1 THEN b0 = %True ELSE RESET b0 IF (ISFALSE (b MOD 2)) OR (ISFALSE (b MOD 3)) THEN b0 = %True END IF IF c = 1 THEN c0 = %True ELSE RESET c0 IF (ISFALSE (c MOD 2)) OR (ISFALSE (c MOD 3)) THEN c0 = %True END IF IF a0 + b0 + c0 = 3 THEN ' in de veralgemeende formulering: sum = (RiHi(0).coefa * a ) + (Rihi(0).coefb * b) + (RiHi(0).coefc * c) ' we kunnen een voorwaarde invoeren om de aantallen te beperken en de indeling in groepen ' te vereenvoudigen: IF (ISFALSE sum MOD 2) OR (ISFALSE sum MOD 3) THEN IF ISFALSE sum MOD FIX(LOG2(sum)-2) THEN ' 02.04.2012 INCR RiHi(sum).aantal ' bevat het aantal maal dat deze ticksum voorkomt IF RiHi(sum).aantal > MaxPermuts THEN MaxPermuts = RiHi(sum).aantal IF sum > RiHiSize THEN RiHiSize = sum RiHi(sum).Tiks(RiHi(sum).aantal-1).a = a RiHi(sum).Tiks(RiHi(sum).aantal-1).b = b RiHi(sum).Tiks(RiHi(sum).aantal-1).c = c END IF END IF END IF END IF NEXT c NEXT b NEXT a ' cleanup: REDIM PRESERVE RiHi(RiHiSize) AS STATIC RiHi_Type ' for debug: IF logflag THEN logfile "possible ticksums and number of occurences for 3-value sets" logfile "Minimum size for number of permutations: " & STR$(MaxPermuts) logfile "Minimumsize for RiHi(): " & STR$(RiHiSize) END IF FOR i = 0 TO UBOUND(RiHi) RiHi(i).coefa = RiHi(0).coefa RiHi(i).coefb = RiHi(0).coefb RiHi(i).coefc = RiHi(0).coefc IF RiHi(i).aantal THEN IF ISFALSE MinTiks THEN MinTiks = i IF logflag THEN logfile STR$(i) & " occurs " & STR$(RiHi(i).aantal) & " times" FOR cnt = 0 TO RiHi(i).aantal -1 logfile " " & STR$(RiHi(i).tiks(cnt).a) & STR$(RiHi(i).tiks(cnt).b) & STR$(RiHi(i).tiks(cnt).c) NEXT cnt END IF END IF NEXT i tog = %True IF logflag THEN logfile "Function 3 initialized " & FUNCNAME$ logfile "Mimimum value for number of ticks : " & STR$(MinTiks) END IF END IF ' normal runtime mode after initialisation: SELECT CASE RiHi(nrtiks).aantal 'sums(nrtiks) CASE %False ' check whether the requested number of ticks has a solution for 1/f FUNCTION = VARPTR(RiHi(nrtiks).Tiks(0)) EXIT FUNCTION CASE 1 ' here we should return the only one answer... FUNCTION = VARPTR(RiHi(nrtiks).Tiks(0)) CASE ELSE ' here filters can be applied as well. DO cnt = (RND(1) ^ SQR(2)) * (RiHi(nrtiks).aantal-1) ' 1/f random distribution LOOP UNTIL RiHi(nrtiks).Tiks(cnt).a > 0 FUNCTION = VARPTR(RiHi(nrtiks).Tiks(cnt)) END SELECT END FUNCTION FUNCTION FindBarUnits4 (BYVAL nrtiks AS LONG) EXPORT AS DWORD ' returns a pointer to a TwelveWords_Type ' find positive integer solutions for the equation '(12a + 6b + 4c + 3d) MOD grouping = %False for the minimum value of nrtiks. ' this maps the 1/f distribution on a set of only four elements. ' The histogram is perfect with these coefficients. ' note that all variables must have a different value for valid solutions. ' 29.03.2012: debug o.k. ' to do: sort abcd groups in inverse order of likelyness ' here we have overflows in the a,b,c,d arrays. There are up to 113 permutations for some sums... ' we need sharper filters. LOCAL a,b,c,d, sum, i , a0, b0,c0, d0 AS LONG STATIC tog, cnt AS INTEGER STATIC dif AS SINGLE ' bepaalt het minimaal verschil geeist tussen de coefficienten STATIC RiHisize, MaxPermuts, MinTiks AS DWORD LOCAL pTw AS TwelveWord_Type PTR 'logfile FUNCNAME$ + STR$(nrtiks) IF ISFALSE tog THEN ' logfile " " + funcname$ + " reinit" ' in this case we calculate the entire lookup DIM RiHi(375) AS STATIC RiHi_Type ' size 374 is minimum required dif = 0.25 ' 25% verschil tussen de waarden pTw = Get1FCoefs (4) ' in g_indep.dll RiHi(0).coefa = @pTw.a '12 RiHi(0).coefb = @pTw.b '6 RiHi(0).coefc = @pTw.c '4 RiHi(0).coefd = @pTw.d '3 RESET cnt, RiHiSize FOR a = 1 TO 16 '32 van een zestiende tot een hele noot, dit zijn het aantal tiks waaruit de waarde van a bestaat. FOR b = 1 TO 16 '32 FOR c = 1 TO 16 FOR d = 1 TO 16 ' IF a <> b AND b <> c AND a <> c AND a <> d AND b<> d AND c <> d THEN ' veralgemeender geformuleerd. IF (ABS(a-b) > MAX(a,b)* dif) AND (ABS(a-c) > MAX(a,c) * dif) AND (ABS(a-d) > MAX(a,d) * dif) AND (ABS(b-c) > MAX(b,c) * dif) AND (ABS(b-d) > MAX(b,d) * dif) AND (ABS(c-d) > MAX(c,d) * dif) THEN ' more conditions required: ' a,b,c,d if different from 1, should be divisible by 2 or 3 IF a = 1 THEN a0 = %True ELSE RESET a0 IF (ISFALSE (a MOD 2)) OR (ISFALSE (a MOD 3)) THEN a0 = %True END IF IF b = 1 THEN b0 = %True ELSE RESET b0 IF (ISFALSE (b MOD 2)) OR (ISFALSE (b MOD 3)) THEN b0 = %True END IF IF c = 1 THEN c0 = %True ELSE RESET c0 IF (ISFALSE (c MOD 2)) OR (ISFALSE (c MOD 3)) THEN c0 = %True END IF IF d = 1 THEN d0 = %True ELSE RESET d0 IF (ISFALSE (d MOD 2)) OR (ISFALSE (d MOD 3)) THEN d0 = %True END IF ' een heel scherpe extra konditie zou kunnen zijn te eisen dat: 'IF (a = 2 * b) OR (a = b / 2) OR (a = 3 * b) OR (a = b/3) THEN IF a0 + b0 + c0 + d0 = 4 THEN 'if (abs ((c + d) /2) - d) < 3 then ' strenge te testen konditie 'IF ABS(c-d) > 3 THEN 'sum = (2 * a) + b , of in de veralgemeende formulering: sum = (RiHi(0).coefa * a ) + (Rihi(0).coefb * b) + (RiHi(0).coefc * c) + (RiHi(0).coefd * d) ' we kunnen een voorwaarde invoeren om de aantallen te beperken en de indeling in groepen ' te vereenvoudigen: IF (ISFALSE sum MOD 2) OR (ISFALSE sum MOD 3) THEN IF ISFALSE sum MOD FIX(LOG2(sum)-2) THEN ' 02.04.2012 INCR RiHi(sum).aantal ' bevat het aantal maal dat deze ticksum voorkomt IF RiHi(sum).aantal > MaxPermuts THEN MaxPermuts = RiHi(sum).aantal IF sum > RiHiSize THEN RiHiSize = sum RiHi(sum).Tiks(RiHi(sum).aantal-1).a = a RiHi(sum).Tiks(RiHi(sum).aantal-1).b = b RiHi(sum).Tiks(RiHi(sum).aantal-1).c = c RiHi(sum).Tiks(RiHi(sum).aantal-1).d = d END IF END IF END IF 'END IF END IF NEXT d NEXT c NEXT b NEXT a ' cleanup: REDIM PRESERVE RiHi(RiHiSize) AS STATIC RiHi_Type ' for debug: IF logflag THEN logfile "possible ticksums and number of occurences for 4-value sets" logfile "Minimum size for nr.of permutations: " & STR$(MaxPermuts) logfile "Minimum sizing for RiHi(): " & STR$(RiHiSize) END IF FOR i = 0 TO UBOUND(RiHi) RiHi(i).coefa = RiHi(0).coefa RiHi(i).coefb = RiHi(0).coefb RiHi(i).coefc = RiHi(0).coefc RiHi(i).coefd = RiHi(0).coefd IF RiHi(i).aantal THEN IF ISFALSE MinTiks THEN Mintiks = i IF logflag THEN logfile STR$(i) & " occurs " & STR$(RiHi(i).aantal) & " times" FOR cnt = 0 TO RiHi(i).aantal -1 logfile " " & STR$(RiHi(i).tiks(cnt).a) & STR$(RiHi(i).tiks(cnt).b) & STR$(RiHi(i).tiks(cnt).c) & STR$(RiHi(i).tiks(cnt).d) NEXT cnt END IF END IF NEXT i tog = %True IF logflag THEN logfile "Function 4 initialized " & FUNCNAME$ logfile "Mimimum value for number of ticks : " & STR$(MinTiks) END IF END IF ' normal runtime mode after initialisation: ' logfile "aantal" + str$(RiHi(nrtiks).aantal) + " @" + funcname$ SELECT CASE RiHi(nrtiks).aantal CASE %False ' check whether the requested number of ticks has a solution for 1/f FUNCTION = VARPTR(RiHi(nrtiks).Tiks(0)) EXIT FUNCTION CASE 1 ' here we should return the only one answer... FUNCTION = VARPTR(RiHi(nrtiks).Tiks(0)) CASE ELSE ' here filters can be applied as well. DO cnt = (RND(1) ^ SQR(2)) * (RiHi(nrtiks).aantal-1) ' 1/f random distribution LOOP UNTIL RiHi(nrtiks).Tiks(cnt).a > 0 FUNCTION = VARPTR(RiHi(nrtiks).Tiks(cnt)) END SELECT END FUNCTION FUNCTION FindBarUnits5 (BYVAL nrtiks AS LONG) EXPORT AS DWORD ' returns a pointer to a TwelveWords_Type ' find positive integer solutions for the equation '(60a + 30b + 20c + 15d + 12e) MOD grouping = %False for the minimum value of nrtiks. ' this maps the 1/f distribution on a set of only four elements. ' The histogram is perfect with these coefficients. ' note that all variables must have a different value for valid solutions. LOCAL a,b,c,d,e, sum, i , a0, b0,c0, d0, e0 AS LONG STATIC tog, cnt AS INTEGER STATIC dif AS SINGLE ' bepaalt het minimaal verschil geeist tussen de coefficienten STATIC RiHisize, MaxPermuts, MinTiks AS DWORD LOCAL pTw AS TwelveWord_Type PTR IF ISFALSE tog THEN ' in this case we calculate the entire lookup DIM RiHi(1500) AS STATIC RiHi_Type ' size 1464 is minimum required dif = 0.33 '0.25 ' 25% verschil tussen de waarden ' RiHi(0).coefa = 60 ' RiHi(0).coefb = 30 ' RiHi(0).coefc = 20 ' RiHi(0).coefd = 15 ' RiHi(0).coefe = 12 pTw = Get1FCoefs (5) ' in g_indep.dll RiHi(0).coefa = @pTw.a ' = 60 RiHi(0).coefb = @pTw.b '30 RiHi(0).coefc = @pTw.c '20 RiHi(0).coefd = @pTw.d '15 RiHi(0).coefe = @pTw.e '12 RESET cnt, RiHiSize FOR a = 1 TO 16 '32 van een zestiende tot een hele noot, dit zijn het aantal tiks waaruit de waarde van a bestaat. FOR b = 1 TO 16 '32 FOR c = 1 TO 16 FOR d = 1 TO 16 FOR e = 1 TO 16 ' veralgemeend geformuleerd. IF (ABS(a-b) > MAX(a,b)* dif) AND (ABS(a-c) > MAX(a,c) * dif) AND (ABS(a-d) > MAX(a,d) * dif) AND (ABS(b-c) > MAX(b,c) * dif) AND (ABS(b-d) > MAX(b,d) * dif) AND (ABS(c-d) > MAX(c,d) * dif) THEN IF (ABS(a-e) > MAX(a,e) * dif) AND (ABS(b-e) > MAX(b,e) * dif) AND (ABS(c-e) > MAX(c,e) * dif) AND (ABS(d-e) > MAX(d,e) * dif) THEN ' more conditions required: ' a,b,c,d,e if different from 1, should be divisible by 2 or 3 IF a = 1 THEN a0 = %True ELSE RESET a0 IF (ISFALSE (a MOD 2)) OR (ISFALSE (a MOD 3)) THEN a0 = %True END IF IF b = 1 THEN b0 = %True ELSE RESET b0 IF (ISFALSE (b MOD 2)) OR (ISFALSE (b MOD 3)) THEN b0 = %True END IF IF c = 1 THEN c0 = %True ELSE RESET c0 IF (ISFALSE (c MOD 2)) OR (ISFALSE (c MOD 3)) THEN c0 = %True END IF IF d = 1 THEN d0 = %True ELSE RESET d0 IF (ISFALSE (d MOD 2)) OR (ISFALSE (d MOD 3)) THEN d0 = %True END IF IF e = 1 THEN e0 = %True ELSE RESET e0 IF (ISFALSE (e MOD 2)) OR (ISFALSE (e MOD 3)) THEN e0 = %True END IF IF a0 + b0 + c0 + d0 + e0 = 5 THEN ' in de veralgemeende formulering: sum = (RiHi(0).coefa * a ) + (Rihi(0).coefb * b) + (RiHi(0).coefc * c) + (RiHi(0).coefd * d) + (RiHi(0).coefe * e) ' we kunnen een voorwaarde invoeren om de aantallen te beperken en de indeling in groepen ' te vereenvoudigen: IF (ISFALSE sum MOD 2) OR (ISFALSE sum MOD 3) THEN IF ISFALSE sum MOD FIX(LOG2(sum)-2) THEN ' 02.04.2012 INCR RiHi(sum).aantal ' bevat het aantal maal dat deze ticksum voorkomt IF RiHi(sum).aantal > MaxPermuts THEN MaxPermuts = RiHi(sum).aantal IF sum > RiHiSize THEN RiHiSize = sum RiHi(sum).Tiks(RiHi(sum).aantal-1).a = a RiHi(sum).Tiks(RiHi(sum).aantal-1).b = b RiHi(sum).Tiks(RiHi(sum).aantal-1).c = c RiHi(sum).Tiks(RiHi(sum).aantal-1).d = d RiHi(sum).Tiks(RiHi(sum).aantal-1).e = e END IF END IF END IF END IF END IF NEXT e NEXT d NEXT c NEXT b NEXT a ' cleanup: REDIM PRESERVE RiHi(RiHiSize) AS STATIC RiHi_Type ' for debug: IF logflag THEN logfile "possible ticksums and number of occurences for 5-value sets" logfile "Minimum size for nr.of permutations: " & STR$(MaxPermuts) logfile "Minimum sizing for RiHi(): " & STR$(RiHiSize) END IF FOR i = 0 TO UBOUND(RiHi) RiHi(i).coefa = RiHi(0).coefa RiHi(i).coefb = RiHi(0).coefb RiHi(i).coefc = RiHi(0).coefc RiHi(i).coefd = RiHi(0).coefd RiHi(i).coefe = RiHi(0).coefe IF RiHi(i).aantal THEN IF ISFALSE MinTiks THEN MinTiks = i IF logflag THEN logfile STR$(i) & " occurs " & STR$(RiHi(i).aantal) & " times" FOR cnt = 0 TO RiHi(i).aantal -1 logfile " " & STR$(RiHi(i).tiks(cnt).a) & STR$(RiHi(i).tiks(cnt).b) & STR$(RiHi(i).tiks(cnt).c) & STR$(RiHi(i).tiks(cnt).d) & STR$(RiHi(i).tiks(cnt).e) NEXT cnt END IF END IF NEXT i tog = %True IF logflag THEN logfile "Function 5 initialized " & FUNCNAME$ logfile "Mimimum value for number of ticks : " & STR$(MinTiks) END IF END IF ' normal runtime mode after initialisation: SELECT CASE RiHi(nrtiks).aantal CASE %False ' check whether the requested number of ticks has a solution for 1/f FUNCTION = VARPTR(RiHi(nrtiks).Tiks(0)) EXIT FUNCTION CASE 1 ' here we should return the only one answer... FUNCTION = VARPTR(RiHi(nrtiks).Tiks(0)) CASE ELSE ' here filters can be applied as well. DO cnt = (RND(1) ^ SQR(2)) * (RiHi(nrtiks).aantal-1) ' 1/f random distribution LOOP UNTIL RiHi(nrtiks).Tiks(cnt).a > 0 FUNCTION = VARPTR(RiHi(nrtiks).Tiks(cnt)) END SELECT END FUNCTION FUNCTION FindBarUnits6 (BYVAL nrtiks AS LONG) EXPORT AS DWORD ' returns a pointer to a TwelveWords_Type ' find positive integer solutions for the equation '(60a + 30b + 20c + 15d + 12e + 10f) MOD grouping = %False for the minimum value of nrtiks. ' this maps the 1/f distribution on a set of only six elements. ' The histogram is perfect with these coefficients. ' note that all variables must have a different value for valid solutions. ' 30.03.2012: debug gwr. ok. LOCAL a,b,c,d,e,f, sum, i , a0, b0,c0, d0, e0, f0 AS LONG STATIC tog, cnt AS INTEGER STATIC dif AS SINGLE ' bepaalt het minimaal verschil geeist tussen de coefficienten STATIC RiHisize, MaxPermuts, Mintiks AS DWORD LOCAL pTw AS Twelveword_Type PTR ' for coefficient lookup IF ISFALSE tog THEN ' in this case we calculate the entire lookup DIM RiHi(1600) AS STATIC RiHi_Type ' size 1474 is minimum required dif = 0.33 '0.25 ' 33% verschil tussen de waarden pTw = Get1FCoefs (6) ' in g_indep.dll RiHi(0).coefa = @pTw.a ' = 60 RiHi(0).coefb = @pTw.b '30 RiHi(0).coefc = @pTw.c '20 RiHi(0).coefd = @pTw.d '15 RiHi(0).coefe = @pTw.e '12 Rihi(0).coeff = @pTw.f '10 RESET cnt, RiHiSize FOR a = 1 TO 16 '32 van een zestiende tot een hele noot, dit zijn het aantal tiks waaruit de waarde van a bestaat. FOR b = 1 TO 16 '32 FOR c = 1 TO 16 FOR d = 1 TO 16 FOR e = 1 TO 16 FOR f = 1 TO 16 ' veralgemeend geformuleerd. IF (ABS(a-b) > MAX(a,b)* dif) AND (ABS(a-c) > MAX(a,c) * dif) AND (ABS(a-d) > MAX(a,d) * dif) AND (ABS(b-c) > MAX(b,c) * dif) AND (ABS(b-d) > MAX(b,d) * dif) AND (ABS(c-d) > MAX(c,d) * dif) THEN IF (ABS(a-e) > MAX(a,e) * dif) AND (ABS(b-e) > MAX(b,e) * dif) AND (ABS(c-e) > MAX(c,e) * dif) AND (ABS(d-e) > MAX(d,e) * dif) THEN IF (ABS(a-f) > MAX(a,f) * dif) AND (ABS(b-f) > MAX(b,f) * dif) AND (ABS(c-f) > MAX(c,f) * dif) AND (ABS(d-f) > MAX(d,f) * dif) AND (ABS(e-f) > MAX(e,f) * dif) THEN ' more conditions required: ' a,b,c,d,e if different from 1, should be divisible by 2 or 3 IF a = 1 THEN a0 = %True ELSE RESET a0 IF (ISFALSE (a MOD 2)) OR (ISFALSE (a MOD 3)) THEN a0 = %True END IF IF b = 1 THEN b0 = %True ELSE RESET b0 IF (ISFALSE (b MOD 2)) OR (ISFALSE (b MOD 3)) THEN b0 = %True END IF IF c = 1 THEN c0 = %True ELSE RESET c0 IF (ISFALSE (c MOD 2)) OR (ISFALSE (c MOD 3)) THEN c0 = %True END IF IF d = 1 THEN d0 = %True ELSE RESET d0 IF (ISFALSE (d MOD 2)) OR (ISFALSE (d MOD 3)) THEN d0 = %True END IF IF e = 1 THEN e0 = %True ELSE RESET e0 IF (ISFALSE (e MOD 2)) OR (ISFALSE (e MOD 3)) THEN e0 = %True END IF IF f = 1 THEN f0 = %True ELSE RESET f0 IF (ISFALSE (f MOD 2)) OR (ISFALSE (f MOD 3)) THEN f0 = %True END IF IF a0 + b0 + c0 + d0 + e0 + f0 = 6 THEN ' in de veralgemeende formulering: sum = (RiHi(0).coefa * a ) + (Rihi(0).coefb * b) + (RiHi(0).coefc * c) + (RiHi(0).coefd * d) + (RiHi(0).coefe * e) + (RiHi(0).coeff * f) ' we kunnen een voorwaarde invoeren om de aantallen te beperken en de indeling in groepen ' te vereenvoudigen: IF (ISFALSE sum MOD 2) OR (ISFALSE sum MOD 3) THEN IF ISFALSE sum MOD FIX(LOG2(sum)-2) THEN ' 02.04.2012 INCR RiHi(sum).aantal ' bevat het aantal maal dat deze ticksum voorkomt IF RiHi(sum).aantal > MaxPermuts THEN MaxPermuts = RiHi(sum).aantal IF sum > RiHiSize THEN RiHiSize = sum RiHi(sum).Tiks(RiHi(sum).aantal-1).a = a RiHi(sum).Tiks(RiHi(sum).aantal-1).b = b RiHi(sum).Tiks(RiHi(sum).aantal-1).c = c RiHi(sum).Tiks(RiHi(sum).aantal-1).d = d RiHi(sum).Tiks(RiHi(sum).aantal-1).e = e RiHi(sum).Tiks(RiHi(sum).aantal-1).f = f END IF END IF END IF END IF END IF END IF NEXT f NEXT e NEXT d NEXT c NEXT b NEXT a ' cleanup: REDIM PRESERVE RiHi(RiHiSize) AS STATIC RiHi_Type ' for debug: IF logflag THEN logfile "possible ticksums and number of occurences for 6-value sets" logfile "Minimum size for nr.of permutations: " & STR$(MaxPermuts) logfile "Minimum sizing for RiHi(): " & STR$(RiHiSize) END IF FOR i = 0 TO UBOUND(RiHi) RiHi(i).coefa = RiHi(0).coefa RiHi(i).coefb = RiHi(0).coefb RiHi(i).coefc = RiHi(0).coefc RiHi(i).coefd = RiHi(0).coefd RiHi(i).coefe = RiHi(0).coefe RiHi(i).coeff = RiHi(0).coeff IF RiHi(i).aantal THEN IF ISFALSE MinTiks THEN Mintiks = i ' kleinst mogelijke waarde voor nrtiks IF logflag THEN logfile STR$(i) & " occurs " & STR$(RiHi(i).aantal) & " times" FOR cnt = 0 TO RiHi(i).aantal -1 logfile " " & STR$(RiHi(i).tiks(cnt).a) & STR$(RiHi(i).tiks(cnt).b) & STR$(RiHi(i).tiks(cnt).c) & STR$(RiHi(i).tiks(cnt).d) & STR$(RiHi(i).tiks(cnt).e) & STR$(RiHi(i).tiks(cnt).f) NEXT cnt END IF END IF NEXT i tog = %True IF logflag THEN logfile "Function 6 initialized " & FUNCNAME$ logfile "Mimimum value for number of ticks : " & STR$(MinTiks) ' test kgv fuction: DIM a(5) AS WORD AT pTW 'varptr i = kgv_arw (a()) logfile "kgv_arw returns= " & STR$(i) ' =60 END IF END IF ' normal runtime mode after initialisation: SELECT CASE RiHi(nrtiks).aantal CASE %False ' check whether the requested number of ticks has a solution for 1/f FUNCTION = VARPTR(RiHi(nrtiks).Tiks(0)) EXIT FUNCTION CASE 1 ' here we should return the only one answer... FUNCTION = VARPTR(RiHi(nrtiks).Tiks(0)) CASE ELSE ' here filters can be applied as well. DO cnt = (RND(1) ^ SQR(2)) * (RiHi(nrtiks).aantal-1) ' 1/f random distribution LOOP UNTIL RiHi(nrtiks).Tiks(cnt).a > 0 FUNCTION = VARPTR(RiHi(nrtiks).Tiks(cnt)) END SELECT END FUNCTION FUNCTION ggd (BYVAL a AS DWORD, BYVAL b AS DWORD) EXPORT AS DWORD ' algoritme van Euclides voor de berekening van de grootste gemene deler van twee getallen ' 31.03.2012: gwr - checked o.k. ' in english, this function ought to be called gcd(), greatest common denominator LOCAL rest AS DWORD ' changed to dword 30.09.2018 , was INTEGER IF a < b THEN SWAP a,b ' a must be > b IF ISFALSE b THEN FUNCTION = %False : EXIT FUNCTION DO rest = a MOD b a = b b = rest LOOP WHILE rest '<> 0 FUNCTION = a END FUNCTION FUNCTION kgv (BYVAL a AS DWORD, BYVAL b AS DWORD) EXPORT AS DWORD ' funktie voor de berekening van het kleinste gemeen veelvoud van twee getallen ' 31.03.2012: gwr - checked o.k. LOCAL g AS DWORD g = ggd (a,b) IF ISFALSE g THEN FUNCTION = %False : EXIT FUNCTION FUNCTION = ABS(a * b) / g END FUNCTION FUNCTION kgv_arw (BYREF ar() AS WORD) EXPORT AS DWORD ' to be checked - ar() should not contain any zero's !!! LOCAL i AS WORD LOCAL g AS DWORD g = kgv (ar(0), ar(1)) FOR i = 2 TO UBOUND(ar) g = kgv (g, ar(i)) NEXT i FUNCTION = g END FUNCTION FUNCTION kgv_ardw (BYREF ar() AS DWORD) EXPORT AS DWORD ' to be checked LOCAL i AS WORD LOCAL g AS DWORD g = kgv (ar(0), ar(1)) FOR i = 2 TO UBOUND(ar) g = kgv (g, ar(i)) NEXT i FUNCTION = g END FUNCTION FUNCTION HistoTest (BYREF His AS HistoGram_Type, BYVAL newidx AS DWORD) EXPORT AS LONG ' "Good Music" 1/f distribution test procedure. ' newidx is the index in dta() for which we have to check the conditions. ' the function returns %True or %False. ' His.fault returns a weighted value for the quality of the distribution ' this is a multi-instance function in g_indep.dll ' For generalisation we created a specific HistoGram_Type that we can pass. If we make arrays of them we can easily have multiple instances. ' Note that it is up to the user to dimension the required arrays as well as to pass the pointers prior to using this function. ' always first set the pointers: LOCAL fault, sum, oldfout, oldfault AS SINGLE LOCAL i,j, minval, maxval, tekort AS LONG IF newidx > His.nr THEN EXIT FUNCTION ' note: newidx must be <= His.nr !!!! IF His.nr < 1 THEN EXIT FUNCTION ' there must be at least 2 elements. IF ISFALSE His.tog THEN ' first check whether valid pointers have been passed: IF His.pDta THEN REDIM Dta(His.nr) AS LONG AT His.pDta ELSE Warning "NUL-pointer passed for dta in " & FUNCNAME$ EXIT FUNCTION END IF IF His.pHisto THEN REDIM Histo(His.nr) AS DWORD AT His.pHisto ' can only have positive values. ' histo contains the number of occurences of any element in dta() ELSE Warning "NUL-pointer passed for Histo in " & FUNCNAME$ EXIT FUNCTION END IF IF His.pFout THEN REDIM Fout(His.nr) AS SINGLE AT His.pFout ' contains the error of each value in the histogram ELSE Warning "NUL-pointer passed for Fout in " & FUNCNAME$ EXIT FUNCTION END IF ' initialize by filling in the Soll data proportions: Histo(0) = (His.nr + 1) * (His.nr + 1) j = 2 FOR i = 1 TO His.nr Histo(i) = Histo(0) / j ' /2, /3, /4, /5, /6,.... ' should be a correct 1/f distribution INCR j NEXT i ' initialize the error array: (using the algorithm) GOSUB Histo_Recalc ' fout() ought to have (near) zero errors on init. His.tog = %True ELSE ' make sure every instance of the function always refers to its own dataset. REDIM Dta(His.nr) AS LONG AT His.pDta REDIM Histo(His.nr) AS DWORD AT His.pHisto REDIM Fout(His.nr) AS SINGLE AT His.pFout END IF GOSUB Histo_Recalc ' we add the new element, and if unacceptable, remove it again. ' fout(newidx) contains the error without the new element ' fault contains the total weighted error of the distribution oldfout = fout(newidx) oldfault = fault ' now add the new element: INCR Histo(newidx) GOSUB Histo_Recalc ' and evaluate the result: IF ABS(fault) =< ABS(oldfault) THEN ' accept result if the global error has decreased his.fault = fault 'total error = fault FUNCTION = %True EXIT FUNCTION END IF SELECT CASE ABS(fout(newidx)) CASE <= ABS(oldfout) ' accept result his.fault = fault 'total error = fault FUNCTION = %True EXIT FUNCTION CASE <= His.tol ' 0.1 = 10% error accepted His.fault = fault ' total error = fault FUNCTION = %True EXIT FUNCTION CASE ELSE ' reject the new element DECR Histo(newidx) ' remove it from histo FUNCTION = %False END SELECT GOSUB Histo_Recalc His.Fault = Fault 'toterr = Fault EXIT FUNCTION Histo_Recalc: ' recalculate the total error-array: 'RESET minval ' was bug ' here we try to detect the note index for which the fault is minimal. ' this note ought to be the most adequate note to be played. However, it will always increase total error. ' but... we could also suggest the note that would contribute most to minimalisation of the error. ' Note however, that we can only add notes, never remove any from the histogram. minval = LargestLNG '&H7FFFFFFF 'reset maxval His.sug = -1 FOR i = 1 TO his.nr fout(i) = Histo(i) - (Histo(0)/(i+1)) ' 1/f distributie test - this means that fout(0) is always false!!! 'IF fout(i) < minval THEN minval = fout(i): His.sug = i - changed to: ' IF ABS(fout(i)) < minval THEN minval = ABS(Fout(i)) : His.sug = i 'IF fout(i) > maxval then maxval = fout(i) : His.sug = i dit zou de fout alleen maar vergroten... 'if fout(i) > maxval then maxval = fout(i) : His.sug = i -1 -ook fout 'iF fout(i) < minval then minval = Fout(i) 'IF fout(i) > maxval then maxval = Fout(i) fout(i) = fout(i) / Histo(i) ' normalize NEXT i ' we recalculate the fault on fout(0) as: fout(0) = Histo(0) - (Histo(1) * 2) ' 1/f rule ' IF ABS(fout(0)) < minval THEN minval = ABS(fout(0)): His.sug = 0 ' ABS added 16.04.2012 fout(0) = fout(0) / Histo(0) ' calculate total degree of error in the distribution as a weighted sum: j = 1 RESET fault, sum FOR i = His.nr TO 0 STEP -1 fault = fault + (fout(i) * j) sum = sum + j INCR j ' 1/f NEXT i ' normalize: fault = fault / sum ' new algo 17.04.2012: find an index value to suggest: IF histo(0) < (histo(1) * 2) THEN his.sug = 0 ELSE ' zoek het grootste tekort in de distributie RESET maxval FOR i = 0 TO his.nr - 1 tekort = (histo(i) * (i/(i+1)) ) - histo(i+1) IF tekort > maxval THEN maxval = tekort : his.sug = i + 1 NEXT i END IF RETURN END FUNCTION SUB Decr0 (i AS BYTE) EXPORT ' decrement with ceiling to 0 IF i > %False THEN DECR i ELSE i = %False END SUB SUB Analyze_Number (BYVAL n AS DWORD) EXPORT ' writes the output to a file '30/01/2013 sommen van 4 priemen herschreven door xof - nu krijgen we veel meer resultaten, maar de berekening kan ook heel lang duren.. ' to do: checken of dit klopt en dan 3, 5, 7, ... aanpassen met zelfde algoritme LOCAL p,q,r,s,t,u,v,w,x,y,z,a, b, sum, cnt, m AS DWORD m = n logfile "Analyse van het getal " & STR$(n), STR$(n)& "_" logfile "Hexadecimaal: " & HEX$(m) logfile "Octaal: " & OCT$(m) logfile "Binair: " & BIN$(m) logfile "Ontbinding in priemfaktoren:" ' ontbinding van een getal in priemfaktoren p = nextprime(1) ' = 2 DO RESET cnt DO IF n MOD p = 0 THEN n = n / p INCR cnt ELSE EXIT DO END IF LOOP 'until n mod p ' output cnt en p naar file IF cnt THEN logfile STR$(cnt) & " x " & STR$(p) & " or " & STR$(p) & "^" & STR$(cnt) END IF p = nextprime(p) LOOP UNTIL p >= n logfile "1 x " & STR$(n) & " or " & STR$(n) & "^ 1" ' ontbinding van een getal als som van priemgetallen logfile "Ontbinding van het getal als som van verschillende priemgetallen:" p = 1 'nextprime(1) q = 2 'nextprime(p) r = nextprime(q) s = nextprime(r) t = nextprime(s) ' sommen van 2 priemen: - 2011+2 wordt niet gevonden hier... IF m MOD 2 = 0 THEN ' logfile "Sums of 2 prime numbers:" DO DO sum = p + q IF sum = m THEN logfile STR$(p) & " + " & STR$(q) & " =" & STR$(m) EXIT DO END IF q = nextprime(q) LOOP UNTIL sum >= m p = nextprime(p) q = nextprime(p) LOOP UNTIL p + q > m END IF ' sommen van 3 priemen 'xof remark: ontbindingen waarin twee keer hetzelfde getal voorkomt zijn ongewild? (bv 11 = 3 + 3 + 5) zal op deze manier nietr gevonden worden.. ' ook verschillende oplossingen met eenzelfde begingetal worden geskipt - of zijn die er niet?? logfile "Sums of 3 prime numbers:" p = 1 'nextprime(1) q = 2 'nextprime(p) r = nextprime(q) s = nextprime(r) t = nextprime(s) DO DO DO sum = p+q+r ' logfile str$(p) + str$(q) + str$(r) + "-" +str$(p+q+r) IF sum = m THEN logfile STR$(p) & " + " & STR$(q) & " + " & STR$(r) & " =" & STR$(m) EXIT DO END IF q = nextprime(q) LOOP UNTIL sum > m r = nextprime(r) q = nextprime(p) 'added by xof 30/01/2013 LOOP UNTIL sum >= m p = nextprime(p) q = nextprime(p) r = nextprime(q) LOOP UNTIL p+q+r>m ' sommen van 4 priemen ' IF m MOD 2 = 0 THEN 'was bug! oplossingen met 2 werden niet gevonden logfile "Sums of 4 prime numbers:" '' p = 1 'nextprime(1) '' q = 2 'nextprime(p) '' r = nextprime(q) '' s = nextprime(r) '' t = nextprime(s) '' DO '' DO '' DO '' DO '' sum = p+q+r+s '' ' logfile str$(p) + str$(q)+str$(r) + str$(s) + "->" + str$(sum) '' IF sum = m THEN '' logfile STR$(p) & " + " & STR$(q) & " + " & STR$(r) & " + " & STR$(s) & " =" & STR$(m) '' EXIT DO '' END IF '' q = nextprime(q) '' LOOP UNTIL sum > m '' ' logfile "--- next r ---" '' r = nextprime(r) '' q = nextprime(p) 'xof als p hier r wordt, missen we resultaten '' s = nextprime(r) 'xof '' sum = p + q + r + s 'xof '' LOOP UNTIL sum >= m '' ' logfile "--- next q ---" '' p = nextprime(p) '' q = nextprime(p) 'xof '' r = nextprime(q) 'xof '' s = nextprime(r) 'xof '' sum = p + q + r + s '' LOOP UNTIL sum >=m '' ' q = nextprime(p) '' q = nextprime(q) 'xof '' r = nextprime(q) '' s = nextprime(r) '' LOOP UNTIL p+q+r+s>m '' ' END IF 'new try for all solutions but no doubles by xof 'this delivers much more solutions than the version above.. and takes several mintues to compute for m = 2013. p = 1 q = 2 r = nextprime(q) s = nextprime(r) DO DO DO DO sum = p + q + r + s ' logfile str$(p) + str$(q)+str$(r) + str$(s) + "->" + str$(sum) IF sum = m THEN logfile STR$(p) & " + " & STR$(q) & " + " & STR$(r) & " + " & STR$(s) & " =" & STR$(m) EXIT DO END IF s = nextprime(s) sum = p + q + r + s LOOP UNTIL sum > m r = nextprime(r) s = nextprime(r) sum = p + q + r + s LOOP UNTIL sum > m q = nextprime(q) r = nextprime(q) s = nextprime(r) sum = p + q + r + s IF (m MOD 2) = 0 AND p > 2 THEN GOTO fives LOOP UNTIL sum > m p = nextprime(p) q = nextprime(p) r = nextprime(q) s = nextprime(r) sum = p + q + r + s LOOP UNTIL sum > m fives: ' sommem van 5 priemgetallen logfile "Sums of 5 prime numbers:" p = 1 'nextprime(1) q = 2 'nextprime(p) r = nextprime(q) s = nextprime(r) t = nextprime(s) DO DO DO DO DO sum = p+q+r+s+ t IF sum = m THEN logfile STR$(p) & " + " & STR$(q) & " + " & STR$(r) & " + " & STR$(s) & " + " & STR$(t) & " =" & STR$(m) EXIT DO END IF q = nextprime(q) LOOP UNTIL sum > m r = nextprime(r) LOOP UNTIL sum >= m s = nextprime(s) LOOP UNTIL sum >=m p = nextprime(p) LOOP UNTIL sum >=m q = nextprime(p) r = nextprime(q) s = nextprime(r) t = nextprime(s) LOOP UNTIL p+q+r+s+ t>m ' sommen van 6 priemgetallen IF m MOD 2 = 0 THEN logfile "Sums of 6 prime numbers:" p = 1 'nextprime(1) q = 2 'nextprime(p) r = nextprime(q) s = nextprime(r) t = nextprime(s) u = nextprime(t) DO DO DO DO DO DO sum = p+q+r+s+ t + u IF sum = m THEN logfile STR$(p) & " + " & STR$(q) & " + " & STR$(r) & " + " & STR$(s) & " + " & STR$(t) & " + " & STR$(u) & " =" & STR$(m) EXIT DO END IF q = nextprime(q) LOOP UNTIL sum > m r = nextprime(r) LOOP UNTIL sum >= m s = nextprime(s) LOOP UNTIL sum >=m t = nextprime(t) LOOP UNTIL sum >=m p = nextprime(p) LOOP UNTIL sum >=m q = nextprime(p) r = nextprime(q) s = nextprime(r) t = nextprime(s) u = nextprime(t) LOOP UNTIL p+q+r+s+ t+u >m END IF ' sommen van 7 priemgetallen logfile "Sums of 7 prime numbers:" p = 1 'nextprime(1) q = 2 'nextprime(p) r = nextprime(q) s = nextprime(r) t = nextprime(s) u = nextprime(t) v = nextprime(u) DO DO DO DO DO DO DO sum = p+q+r+s+ t + u + v IF sum = m THEN logfile STR$(p) & " + " & STR$(q) & " + " & STR$(r) & " + " & STR$(s) & " + " & STR$(t) & " + " & STR$(u) & " + " & STR$(v) &" =" & STR$(m) EXIT DO END IF q = nextprime(q) LOOP UNTIL sum > m r = nextprime(r) LOOP UNTIL sum >= m s = nextprime(s) LOOP UNTIL sum >=m t = nextprime(t) LOOP UNTIL sum >=m u = nextprime(u) LOOP UNTIL sum >=m p = nextprime(p) LOOP UNTIL sum >=m q = nextprime(p) r = nextprime(q) s = nextprime(r) t = nextprime(s) u = nextprime(t) v = nextprime(u) LOOP UNTIL p+q+r+s+ t+u+v >m ' sommen van 8 priemgetallen IF m MOD 2 = 0 THEN logfile "Sums of 8 prime numbers:" p = 1 'nextprime(1) q = 2 'nextprime(p) r = nextprime(q) s = nextprime(r) t = nextprime(s) u = nextprime(t) v = nextprime(u) w = nextprime(v) DO DO DO DO DO DO DO DO sum = p+q+r+s+t+u+v+w IF sum = m THEN logfile STR$(p) & "+" & STR$(q) & "+" & STR$(r) & "+" & STR$(s) & "+" & STR$(t) & "+" & STR$(u) & "+" & STR$(v) & "+" & STR$(w) & " =" & STR$(m) EXIT DO END IF q = nextprime(q) LOOP UNTIL sum > m r = nextprime(r) LOOP UNTIL sum >= m s = nextprime(s) LOOP UNTIL sum >=m t = nextprime(t) LOOP UNTIL sum >=m u = nextprime(u) LOOP UNTIL sum >=m v = nextprime (v) LOOP UNTIL sum >=m p = nextprime(p) LOOP UNTIL sum >=m q = nextprime(p) r = nextprime(q) s = nextprime(r) t = nextprime(s) u = nextprime(t) v = nextprime(u) w = nextprime(v) LOOP UNTIL p+q+r+s+t+u+v+w >m END IF ' sommen van 9 priemgetallen logfile "Sums of 9 prime numbers:" p = 1 'nextprime(1) q = 2 'nextprime(p) r = nextprime(q) s = nextprime(r) t = nextprime(s) u = nextprime(t) v = nextprime(u) w = nextprime(v) x = nextprime(w) DO DO DO DO DO DO DO DO DO sum = p+q+r+s+t+u+v+w+x IF sum = m THEN logfile STR$(p) & "+" & STR$(q) & "+" & STR$(r) & "+" & STR$(s) & "+" & STR$(t) & "+" & STR$(u) & "+" & STR$(v) & "+" & STR$(w) & "+" & STR$(x) &" =" & STR$(m) EXIT DO END IF q = nextprime(q) LOOP UNTIL sum > m r = nextprime(r) LOOP UNTIL sum >= m s = nextprime(s) LOOP UNTIL sum >=m t = nextprime(t) LOOP UNTIL sum >=m u = nextprime(u) LOOP UNTIL sum >=m v = nextprime (v) LOOP UNTIL sum >=m w = nextprime (w) LOOP UNTIL sum >=m p = nextprime(p) LOOP UNTIL sum >=m q = nextprime(p) r = nextprime(q) s = nextprime(r) t = nextprime(s) u = nextprime(t) v = nextprime(u) w = nextprime(v) x = nextprime(w) LOOP UNTIL p+q+r+s+t+u+v+w+x >m ' sommen van 10 priemgetallen IF m MOD 2 = 0 THEN logfile "Sums of 10 prime numbers:" p = 1 'nextprime(1) q = 2 'nextprime(p) r = nextprime(q) s = nextprime(r) t = nextprime(s) u = nextprime(t) v = nextprime(u) w = nextprime(v) x = nextprime(w) y = nextprime(x) DO DO DO DO DO DO DO DO DO DO sum = p+q+r+s+t+u+v+w+x+y IF sum = m THEN logfile STR$(p) & "+" & STR$(q) & "+" & STR$(r) & "+" & STR$(s) & "+" & STR$(t) & "+" & STR$(u) & "+" & STR$(v) & "+" & STR$(w) & "+" & STR$(x) & "+" & STR$(y) & " =" & STR$(m) EXIT DO END IF q = nextprime(q) LOOP UNTIL sum > m r = nextprime(r) LOOP UNTIL sum >= m s = nextprime(s) LOOP UNTIL sum >=m t = nextprime(t) LOOP UNTIL sum >=m u = nextprime(u) LOOP UNTIL sum >=m v = nextprime (v) LOOP UNTIL sum >=m w = nextprime (w) LOOP UNTIL sum >=m x = nextprime(x) LOOP UNTIL sum >= m p = nextprime(p) LOOP UNTIL sum >=m q = nextprime(p) r = nextprime(q) s = nextprime(r) t = nextprime(s) u = nextprime(t) v = nextprime(u) w = nextprime(v) x = nextprime(w) y = nextprime(x) LOOP UNTIL p+q+r+s+t+u+v+w+x+y >m END IF ' sommen van 11 priemgetallen logfile "Sums of 11 prime numbers:" p = 1 'nextprime(1) q = 2 'nextprime(p) r = nextprime(q) s = nextprime(r) t = nextprime(s) u = nextprime(t) v = nextprime(u) w = nextprime(v) x = nextprime(w) y = nextprime(x) z = nextprime(y) DO DO DO DO DO DO DO DO DO DO DO sum = p+q+r+s+t+u+v+w+x+y+z IF sum = m THEN logfile STR$(p) & "+" & STR$(q) & "+" & STR$(r) & "+" & STR$(s) & "+" & STR$(t) & "+" & STR$(u) & "+" & STR$(v) & "+" & STR$(w) & "+" & STR$(x) & "+" & STR$(y)&"+"& STR$(z) & " =" & STR$(m) EXIT DO END IF q = nextprime(q) LOOP UNTIL sum > m r = nextprime(r) LOOP UNTIL sum >= m s = nextprime(s) LOOP UNTIL sum >=m t = nextprime(t) LOOP UNTIL sum >=m u = nextprime(u) LOOP UNTIL sum >=m v = nextprime (v) LOOP UNTIL sum >=m w = nextprime (w) LOOP UNTIL sum >=m x = nextprime(x) LOOP UNTIL sum >= m y = nextprime(y) LOOP UNTIL sum >= m p = nextprime(p) LOOP UNTIL sum >=m q = nextprime(p) r = nextprime(q) s = nextprime(r) t = nextprime(s) u = nextprime(t) v = nextprime(u) w = nextprime(v) x = nextprime(w) y = nextprime(x) z = nextprime(y) LOOP UNTIL p+q+r+s+t+u+v+w+x+y+z >m ' sommen van 12 priemgetallen IF ISFALSE m MOD 2 THEN logfile "Sums of 12 prime numbers:" p = 1 'nextprime(1) q = 2 'nextprime(p) r = nextprime(q) s = nextprime(r) t = nextprime(s) u = nextprime(t) v = nextprime(u) w = nextprime(v) x = nextprime(w) y = nextprime(x) z = nextprime(y) a = nextprime(z) DO DO DO DO DO DO DO DO DO DO DO DO sum = p+q+r+s+t+u+v+w+x+y+z+a IF sum = m THEN logfile STR$(p) & "+" & STR$(q) & "+" & STR$(r) & "+" & STR$(s) & "+" & STR$(t) & "+" & STR$(u) & "+" & _ STR$(v) & "+" & STR$(w) & "+" & STR$(x) & "+" & STR$(y)&"+"& STR$(z) & "+" & STR$(a)&" =" & STR$(m) EXIT DO END IF q = nextprime(q) LOOP UNTIL sum > m r = nextprime(r) LOOP UNTIL sum >= m s = nextprime(s) LOOP UNTIL sum >=m t = nextprime(t) LOOP UNTIL sum >=m u = nextprime(u) LOOP UNTIL sum >=m v = nextprime (v) LOOP UNTIL sum >=m w = nextprime (w) LOOP UNTIL sum >=m x = nextprime(x) LOOP UNTIL sum >= m y = nextprime(y) LOOP UNTIL sum >= m z = nextprime(z) LOOP UNTIL sum >= m p = nextprime(p) LOOP UNTIL sum >=m q = nextprime(p) r = nextprime(q) s = nextprime(r) t = nextprime(s) u = nextprime(t) v = nextprime(u) w = nextprime(v) x = nextprime(w) y = nextprime(x) z = nextprime(y) a = nextprime(z) LOOP UNTIL p+q+r+s+t+u+v+w+x+y+z+a >m END IF ' sommen van 13 priemgetallen logfile "Sums of 13 prime numbers:" p = 1 'nextprime(1) q = 2 'nextprime(p) r = 3 'nextprime(q) s = 5 'nextprime(r) t = 7 'nextprime(s) u = 11 'nextprime(t) v = 13 'nextprime(u) w = 17 'nextprime(v) x = 19 'nextprime(w) y = 23 'nextprime(x) z = 29 'nextprime(y) a = 31 'nextprime(z) b = 37 'nextprime(a) DO DO DO DO DO DO DO DO DO DO DO DO DO sum = p+q+r+s+t+u+v+w+x+y+z+a+b IF sum = m THEN logfile STR$(p)&"+"& STR$(q)& "+" & STR$(r) & "+" & STR$(s) & "+" & STR$(t) & "+" & STR$(u) & "+"_ & STR$(v) & "+" & STR$(w) & "+" & STR$(x) & "+" & STR$(y)&"+"& STR$(z) & "+" & STR$(a)& "+"& STR$(b)&" =" & STR$(m) EXIT DO END IF q = nextprime(q) LOOP UNTIL sum > m r = nextprime(r) LOOP UNTIL sum >= m s = nextprime(s) LOOP UNTIL sum >=m t = nextprime(t) LOOP UNTIL sum >=m u = nextprime(u) LOOP UNTIL sum >=m v = nextprime (v) LOOP UNTIL sum >=m w = nextprime (w) LOOP UNTIL sum >=m x = nextprime(x) LOOP UNTIL sum >= m y = nextprime(y) LOOP UNTIL sum >= m z = nextprime(z) LOOP UNTIL sum >=m a = nextprime(a) LOOP UNTIL sum >= m p = nextprime(p) LOOP UNTIL sum >=m q = nextprime(p) r = nextprime(q) s = nextprime(r) t = nextprime(s) u = nextprime(t) v = nextprime(u) w = nextprime(v) x = nextprime(w) y = nextprime(x) z = nextprime(y) a = nextprime(z) b = nextprime(a) LOOP UNTIL p+q+r+s+t+u+v+w+x+y+z+a+b >m END SUB FUNCTION Chaos (BYVAL k AS SINGLE,BYVAL x AS SINGLE) EXPORT AS SINGLE '30.12.2015: added by gwr. ' for full documentation cfr. Wikipedia under Chaos formula FUNCTION = k * x * ( 1! - x) END FUNCTION FUNCTION Datim () EXPORT AS STRING ' 14.03.2017 added by gwr. FUNCTION = RIGHT$(DATE$,4) & LEFT$(DATE$,2) & MID$(DATE$,4,2) & LEFT$(TIME$,2) & MID$(TIME$,4,2) & RIGHT$(TIME$,2) END FUNCTION FUNCTION PID (BYVAL sollvalue AS SINGLE, BYVAL seinvalue AS SINGLE, OPT BYVAL kp AS SINGLE, OPT BYVAL ki AS SINGLE, OPT BYVAL kd AS SINGLE) EXPORT AS SINGLE ' 30.05.2017: Generic PID regulation algorithm. - code added by gwr. ' The machine constants have to be passed on the first call only. ' seinvalue is the measured reality value, generaly derived from a sample ' sollvalue is the goal we want to achieve ' The function returns the correction factor for regulation. ' The function should be used in a regulation loop. STATIC propconstant, integrationconstant, differenciationconstant AS SINGLE STATIC oldfout, iterm AS SINGLE LOCAL fout, pterm, dterm AS SINGLE IF kp THEN propconstant = kp IF ki THEN IF ki <> integrationconstant THEN RESET iterm ' reset! integrationconstant = ki END IF IF kd THEN IF kd <> differenciationconstant THEN RESET oldfout ' reset differenciationconstant = kd END IF fout = sollvalue - seinvalue ' calculate the error pterm = propconstant * fout ' proportionality term iterm = iterm + (integrationconstant * fout) ' integration term dterm = differenciationconstant * (fout - oldfout) oldfout = fout FUNCTION = pterm + iterm + dterm ' return value for the PID correction signal END FUNCTION FUNCTION RatNumPat (BYVAL number AS LONG) EXPORT AS STRING LOCAL decimals AS DWORD LOCAL num$, a$, b$ ,pat$ LOCAL teller, deling, rest, n, m, i, j AS LONG ' gwr 29.09.2018 - added to library after debug in PBCC ' debugged o.k. ' number: getal waarvoor we een patroon zoeken in het rationaal getal 1/number ' de funktie retourneert een string met alle cijfers van het zich herhalend patroon ' decimals = number * 2 ' seems we get no trouble doing this. ' no math proof for this however... ' 30.09.2018: There is a proof for the statement that the length of the repeating pattern is always smaller than the number. ' The algorithm used here could be generalised to make a function returning a repeating pattern in any array ' Could be usefull for data mining in audio as well as in gesture recognition. decimals = number * 2 ' first we eliminate a few easy cases to save processing load: SELECT CASE number CASE <3, 4, 5, 8, 10, 100, 1000, 10000 FUNCTION = "" : EXIT FUNCTION ' this also prevents divide by 0 bugs. CASE 3 FUNCTION = "3" : EXIT FUNCTION END SELECT ' perform the division algo step by step to get any number of decimals ' algo - deling , zoals geleerd in de lagere school... teller = 1 deling = teller \ number ' integer divide 'num$ = TRIM$(STR$(deling)) & "," ' voor konventionele notatie ' bij veralgemening kunnen we de komma beter weglaten. num$ = TRIM$(STR$(deling)) DO rest = teller - (deling * number) teller = rest * 10 deling = FIX(teller / number) ' must be 0 to 9 num$ = num$ & TRIM$(STR$(deling)) LOOP UNTIL LEN(num$) > decimals ' this is tested 100% ok 'PRINT num$ ' now we have to find a repetitive pattern in num$ ' there is no mathematical guarantee that this is correct however ' we would have to check that the pattern continues to repeat... ' remove the leading 0 at the start: ' num$ = MID$(num$,2, decimals -1) - not sure about this ' make sure the length is even: 'IF LEN(num$) MOD 2 <> 0 THEN ' num$ = LEFT$(num$,LEN(num$)-1) 'END IF ' solving some easy cases first and fast ' this solves all repeats of one and the same digit. i = ASC("0") j = ASC("9") FOR n = i TO j pat$ = STRING$(number,n) IF INSTR (num$, pat$) THEN IF CHR$(n) = "0" THEN FUNCTION = "" ELSE FUNCTION = CHR$(n) END IF EXIT FUNCTION END IF NEXT n n = 2 ' length 1 is already solved here m = 1 ' starting at the beginning of the string 'xof vraagt zich af: als het herhaalde stuk an de vorm "abcbc" was, zou onderstaande code dan niet "bc" herkennen als herhaalde subsring in plaats van het geheel? DO ' lus met m n = 2 DO ' lus met n a$ = MID$(num$,m,n) b$ = MID$(num$,m+ n,n) IF a$ = b$ THEN FUNCTION = a$ : EXIT FUNCTION INCR n IF n > number THEN EXIT DO ' dit is bewezen... ' we kunnen zelfs versnellen door IF n > Phi(number) te stellen. LOOP INCR m ' startpositie doorschuiven IF m > LEN(num$)/2 THEN EXIT DO LOOP ' if we get here, no pattern was found... FUNCTION = "no pattern... error" END FUNCTION FUNCTION Phi (BYVAL n AS LONG) EXPORT AS LONG ' 30.09.2018 - gwr ' calculation of Eulers totient, using a Fourier transform of the ggd's. ' math source for this: Wikipedia entry on Euler totient. ' Pi2 = ATN(1) * 4 LOCAL k, g, retval AS LONG LOCAL j AS DOUBLE FOR k = 1 TO n g = ggd (k,n) j = j + (g * COS(Pi2 * k / n)) NEXT k IF j - INT(j) > 0 THEN retval = INT(j) + 1 ELSE retval = INT(j) ' afronding naar omhoog FUNCTION = retval END FUNCTION ' functions for conversion of numbers in the duodecimal system FUNCTION DuoDec2Dec (num$) EXPORT AS QUAD ' 01.10.2018 - gwr ' suppose duodec using the digits 0,1,2,3,4,5,6,7,8,9,A,B ' debugged and tested o.k. ' however, this can overflow if VAL(num$) > quad ( 2^64 -1) LOCAL n, m AS QUAD LOCAL i , L AS LONG LOCAL d$ L = LEN(num$) 'ok ' so size is 10^L m = 0 FOR i = L TO 1 STEP -1 ' i and L should be long, not dword!!! d$ = MID$(num$,i,1) SELECT CASE d$ CASE "0" n = 0 CASE "1" TO "9" n = VAL(d$) * (12 ^ (L-i)) CASE "A" n = 10 * (12 ^ (L-i)) CASE "B" n = 11 * (12 ^ (L-i)) CASE ELSE warning d$ & " is not a valid digit in duodecimal" END SELECT m = m + n NEXT i FUNCTION = m END FUNCTION FUNCTION DuoDec (BYVAL n AS QUAD) EXPORT AS STRING ' convert decimal to dodeca (base 12 numbers) ' 01.20.2018 tested and debugged gwr. LOCAL r AS QUAD LOCAL ret$ r = 0 ret$ = "" DO r = n MOD 12 SELECT CASE r CASE 0 TO 9 ret$ = TRIM$(STR$(r)) + ret$ CASE 10 ret$ = "A" + ret$ CASE 11 ret$ = "B" + ret$ END SELECT n = n \ 12 LOOP UNTIL n = 0 FUNCTION = ret$ END FUNCTION FUNCTION DecStr2Duodec (n AS STRING) EXPORT AS STRING ' 01.10.2018 - gwr ' convert decimal-string to dozenal numbers (base 12 numbers) ' to be further debugged ' this should handle almost unlimited numbers. LOCAL q,b,r, L,m AS QUAD LOCAL i, j AS LONG LOCAL ret$ IF n = "" OR n = "0" THEN FUNCTION = "0": EXIT FUNCTION L = LEN(n) IF L < 18 THEN ' no need to handle long strings here... FUNCTION = DuoDec(VAL(n)) EXIT FUNCTION ' debugged o.k. END IF DIM segs(0 TO (L/18) + 1) AS STRING ' limit the chunks to the size quads can handle i = 0 j = 1 DO segs(i) = MID$(n,j,18) ' left to right... j = j + 18 INCR i LOOP UNTIL i > UBOUND(segs) segs(UBOUND(segs)) = TRIM$(segs(UBOUND(segs))) ' now we use the numeric conversion algorithm for quads ret$= "" FOR i = 0 TO UBOUND(segs) q = VAL(segs(i)) ' quad r = q MOD 12 ' rest, over te dragen q = q \ 12 ' integer divide ret$ = ret$ + DuoDec(q) NEXT i IF r > 0 THEN ret$ = ret$ + Duodec(r) ' ? IF LEFT$(ret$,1)= "0" THEN ret$ = MID$(ret$,2,LEN(ret$)-1) FUNCTION = ret$ END FUNCTION '[EOF]