'*************************************************************************** '* Fast Fourier Transforms * '* written In PowerBasic PB32 * '* Dr.Godfried-Willem Raes * '* 1997 * '*************************************************************************** '15.08.1997: first version '18.08.1997: Graphic features added '29.08.1997: revisited '14.09.1997: comments added for teaching $CPU 80386 $DYNAMIC $COMPILE EXE "C:\PB32\FFT\EXE\FFT.EXE" $DEBUG MAP ON $DEBUG PBDEBUG ON $DEBUG PATH ON $DIM ARRAY $FLOAT NPX $OPTIMIZE SPEED 'Variable types (Structures in C) declarations: TYPE Kompleks real as Single imag as Single END TYPE TYPE Polar mag as Single ang as Single END TYPE 'Constants - PB32 allows only integer constants! %True = -1 %False = 0 'Procedure declares: ******************************************************** 'group 1: Fourier Transforms DECLARE SUB Samp2Spek (SampBuf() AS Single, Spectrum() AS Kompleks) ' returns the spectrum as an array of complex numbers ' (real + imaginary) for a bipolar frequency range. ' This is the complete mathematical transform and this ' procedure is the very nucleus of the demonstration. DECLARE SUB Samp2PowSpek (SampBuf() AS Single, PwSp() AS Single) ' returns the power spectrum components as a ' bipolar array PwSp() including negative frequency components DECLARE SUB Sm2Pw (SampBuf() AS Single, PwSp() As Single) :' experiment ' returns the sum of positive and negative frequency ' components in a single unipolar array PwSp() DECLARE SUB Sm2PosPw (SampBuf() AS Single, PwSp() AS Single) ' returns only the values of the positive frequency ' components in a unipolar array PwSp() DECLARE SUB Sm2NegPw (SampBuf() AS Single, PwSp() AS Single) ' returns only the power values of the negative frequency ' components in the spectrum as a unipolar positive array in PwSp DECLARE SUB Sm2Midi (SampBuf() AS Single, MuSp() AS Integer) ' returns a chromatic power distribution transform in ' midi-note increments. ' The power-values have midi range 0-127. ' The strongest spectral component is always 127. 'group 2: Waveformgenerators DECLARE SUB SineWave (SampBuf() AS Single, Freq AS Single) DECLARE SUB SawTooth (SampBuf() AS Single, Freq AS Single) DECLARE SUB SquareWave (SampBuf() AS Single, Freq AS Single) DECLARE SUB TopHat (SampBuf() AS Single, Wdth AS Integer) DECLARE SUB SincWave (SampBuf() AS Single, Freq AS Single) DECLARE SUB NoiseWave (SampBuf() AS Single, PrbExp AS Single) DECLARE SUB GaussWave (SampBuf() AS Single, Wdth AS Single) DECLARE SUB DiracWave (SampBuf() AS Single, Points AS Single) DECLARE SUB ExponWave (SampBuf() AS Single, Wdth AS Single) DECLARE SUB LeesSampleFile (SampBuf() AS Single, kf AS String) DECLARE SUB SchrijfSampleFile (SampBuf() AS Single, kf AS String) 'group 3: Display procedures DECLARE SUB ShowIO (SampBuf() AS Single, PwSp() AS Single) DECLARE SUB ShowWave (SampBuf() AS Single) DECLARE SUB ShowSpectrum (Spectrum() AS Kompleks) DECLARE SUB ShowPwSp (PwSp() AS Single) DECLARE SUB ShowReal (Spectrum() AS Kompleks) DECLARE SUB ShowImag (Spectrum() AS Kompleks) DECLARE SUB ShowRealImag (Spectrum() AS Kompleks) DECLARE SUB ShowMidiSpectrum (SampBuf() AS Single, MuSp() AS Integer) DECLARE FUNCTION Fraster% (i AS Integer) DECLARE SUB BottomLine () DIM FFTSize AS Integer DIM i AS Integer DIM j AS Integer DIM SampBuf(-1 TO 1) AS Single DIM Spectrum(-1 TO 1) AS Kompleks DIM PwSp(-1 TO 1) AS Single DIM MuSp(0 TO 127) AS Integer SCREEN 12 CLS 0 '**************************************************************************** MainCode: SCREEN 12 CLS 0 LOCATE 1,5: PRINT "Fourier transform demonstration program"; LOCATE 2,5: PRINT "---------------------------------------"; LOCATE 3,5: PRINT "by Prof.Dr.Godfried-Willem RAES, 1997 "; LOCATE 7, 10: INPUT "Number of samples (power of 2)?"; FFTSize IF FFTSize=0 THEN FFTSize = 1024 :' default ' FFTSize must be power of 2: 64,128,256,512,1024 ' this is the same as the number of samples LOCATE 10,5: PRINT "Select waveform for input sample:"; LOCATE 11,10: PRINT "1.- Sinewave "; LOCATE 12,10: PRINT "2.- Sawtooth "; LOCATE 13,10: PRINT "3.- Squarewave "; LOCATE 14,10: PRINT "4.- TopHat "; LOCATE 15,10: PRINT "5.- Sin(x)/x "; LOCATE 16,10: PRINT "6.- Noise "; LOCATE 17,10: PRINT "7.- Gauss curve "; LOCATE 18,10: PRINT "8.- Dirac spikes "; LOCATE 19,10: PRINT "9.- Exponential wave "; LOCATE 20,10: PRINT "A.- Read Wave File "; LOCATE 21,10: PRINT "B.- Beta Curve Wave "; ' drawn wave to be added... ' ring-type sample memory buffer to be added LOCATE 22,10: PRINT "ESC.- QUIT" DO k$= INKEY$ IF k$= CHR$(27) THEN END wavekeuze% = VAL("&H"+ k$) LOOP UNTIL wavekeuze% >=1 AND wavekeuze% <=12 SELECT CASE wavekeuze% CASE 1,2,3 CLS 0 :' sine, saw, square REDIM SampBuf(0 TO FFTSize-1) AS Single LOCATE 10,5: PRINT "Frequency Range?"; LOCATE 12,10: PRINT "Minimum Frequency = 1 Hz"; LOCATE 13,10: PRINT "Maximum Frequency ="; FFTSize /2;"Hz"; LOCATE 15,10: INPUT "Start frequency?"; startf! LOCATE 16,10: INPUT "Stop frequency?"; stopf! CASE 4 ' tophat - uses negative time-scale! CLS 0 REDIM SampBuf(-FFTSize\2 TO (FFTSize\2)-1) AS Single LOCATE 10,5: PRINT "Width modulation Range?"; LOCATE 12,10: PRINT "Minimum Width = 2 samples"; LOCATE 13,10: PRINT "Maximum Width ="; FFTSize /2;"samples"; LOCATE 15,10: INPUT "Start width?"; startf! LOCATE 16,10: INPUT "Stop width?"; stopf! CASE 5 ' sin(x)/x - uses negative time-scale CLS 0 REDIM SampBuf(-FFTSize\2 TO (FFTSize\2)-1) AS Single LOCATE 10,5: PRINT "Frequency Range?"; LOCATE 12,10: PRINT "Minimum Frequency = 1 Hz"; LOCATE 13,10: PRINT "Maximum Frequency ="; FFTSize /2;"Hz"; LOCATE 15,10: INPUT "Start frequency?"; startf! LOCATE 16,10: INPUT "Stop frequency?"; stopf! CASE 6 ' noise CLS 0 LOCATE 10,5: PRINT "Parameters ?"; LOCATE 12,10: PRINT "1.- From t0 to tmax - positive time values"; LOCATE 13,10: PRINT "2.- From tmin to tmax - symmetric time values"; LOCATE 14,10: PRINT "Choice ?"; DO noisekeuze$= INKEY$ LOOP UNTIL (noisekeuze$ = "1") OR (noisekeuze$ = "2") noisekeuze% = VAL(noisekeuze$) SELECT CASE noisekeuze% CASE 1 REDIM SampBuf(0 TO FFTSize-1) AS Single CASE 2 REDIM SampBuf((-FFTSize\2) TO (FFTSize\2)-1) AS Single END SELECT LOCATE 16, 10: INPUT "Probability start exponent"; startf! LOCATE 17, 10: INPUT "Probability stop exponent"; stopf! CASE 7 ' Gauss CLS 0 REDIM SampBuf(-FFTSize\2 TO (FFTSize\2)-1) AS Single LOCATE 10,5: PRINT "Width modulation Range?"; LOCATE 12,10: PRINT "Minimum Width =";" 1" ; LOCATE 13,10: PRINT "Maximum Width ="; FFTSize\ 2;" "; LOCATE 15,10: INPUT "Start width?"; startf! LOCATE 16,10: INPUT "Stop width?"; stopf! CASE 8 ' dirac (finite version) CLS 0 REDIM SampBuf(-FFTSize\2 TO (FFTSize\2)-1) AS Single LOCATE 10,5: PRINT "Spike placement?"; LOCATE 12,10: PRINT "Minimum place =";" 0" ; LOCATE 13,10: PRINT "Maximum Width ="; FFTSize\ 2;" "; LOCATE 15,10: INPUT "Start place?"; startf! LOCATE 16,10: INPUT "Stop place?"; stopf! CASE 9 ' exponential wave CLS 0 REDIM SampBuf(-FFTSize\2 TO (FFTSize\2)-1) AS Single LOCATE 10,5: PRINT "Width modulation Range?"; LOCATE 12,10: PRINT "Minimum Slope =";" 1" ; LOCATE 13,10: PRINT "Maximum Slope ="; FFTSize\ 2;" "; LOCATE 15,10: INPUT "Start slope?"; startf! LOCATE 16,10: INPUT "Stop slope?"; stopf! CASE 10 ' read sample data from file on disk CLS 0 LOCATE 10, 5 : PRINT "The data file must be in 16-bit integer format"; LOCATE 11, 5 : PRINT "Each sample, regardless resolution, must occupy"; LOCATE 12, 5 : PRINT "a single two's complement encoded integer"; LOCATE 13, 5: PRINT "The sample data is assumed to be normalized to"; LOCATE 14, 5: PRINT "a lenght of one second"; LOCATE 15, 5: PRINT "The number of samples should be a power of 2"; LOCATE 17, 5: INPUT "File containing sample data ?";FileName$ REDIM SampBuf(0 TO FFTSize - 1) AS Single startf! = 1 stopf!= 1 CASE 11 ' beta functions with 2 parameters REDIM SampBuf(0 TO FFTSize-1) AS Single CASE ELSE END END SELECT '---------------------------------------------------------------------- CLS 0 LOCATE 10, 5: PRINT "Choose Transform type and parameters..."; LOCATE 12, 10: PRINT "1.- Transform to a spectral complex-number array"; LOCATE 13, 10: PRINT "2.- Transform to bipolar frequency power spectrum"; LOCATE 14, 10: PRINT "3.- Transform to summed unipolar frequency spectrum"; LOCATE 15, 10: PRINT "4.- Transform to positive frequency power spectrum"; LOCATE 16, 10: PRINT "5.- Transform to negative frequency power spectrum"; LOCATE 17, 10: PRINT "6.- Transform to chromatic midi-note power spectrum"; LOCATE 18, 10: PRINT "7.- No FFT transform"; LOCATE 19, 10: DO transformkeuze$=INKEY$ IF transformkeuze$= CHR$(27) THEN END LOOP UNTIL transformkeuze$ > "0" AND transformkeuze$ < "8" transformkeuze% = VAL (transformkeuze$) '--------------------------------------------------------------------- CLS 0 LOCATE 10, 5: PRINT "Choose Display mode..."; SELECT CASE transformkeuze% CASE 1 LOCATE 12, 10: PRINT "1.- Show Spectrum from Complex-number array"; LOCATE 13, 10: PRINT "2.- Show Waveform"; LOCATE 16, 10: PRINT "5.- Show real number spectral values only"; LOCATE 17, 10: PRINT "6.- Show imaginary spectral components only"; LOCATE 18, 10: PRINT "7.- Show real&imaginary parts of spectrum"; LOCATE 20, 10: PRINT "9.- Speed calculation of the transform"; LOCATE 21, 10: PRINT "A.- Write Waveform to file"; CASE 2,3,4,5 LOCATE 14, 10: PRINT "3.- Show input waveform and power spectrum"; LOCATE 13, 10: PRINT "2.- Show Waveform"; LOCATE 15, 10: PRINT "4.- Show PowerSpectrum as calculated"; LOCATE 16, 10: PRINT "5.- Show real number spectral values only"; LOCATE 20, 10: PRINT "9.- Speed calculation of the transform"; LOCATE 21, 10: PRINT "A.- Write Waveform to file"; CASE 6 LOCATE 13, 10: PRINT "2.- Show Waveform"; LOCATE 19, 10: PRINT "8.- Show wave + midi-note chromatic transform"; LOCATE 20, 10: PRINT "9.- Speed calculation of the transform"; LOCATE 21, 10: PRINT "A.- Write Waveform to file"; CASE 7 LOCATE 13, 10: PRINT "2.- Show Waveform"; LOCATE 20, 10: PRINT "9.- Speed calculation of the transform"; LOCATE 21, 10: PRINT "A.- Write Waveform to file"; END SELECT DO displaykeuze$=INKEY$ IF displaykeuze$= CHR$(27) THEN END displaykeuze% = VAL ("&H" + displaykeuze$) LOOP UNTIL displaykeuze% > 0 AND displaykeuze% =< 10 starttime! = TIMER count% =0 FOR frq! = startf! TO stopf! STEP 1 SELECT CASE wavekeuze% CASE 1 SineWave SampBuf(), frq! CASE 2 SawTooth SampBuf(), frq! CASE 3 SquareWave SampBuf(), frq! CASE 4 ' here we use the frq! variable for with-modulation TopHat SampBuf(), INT(frq!) CASE 5 SincWave SampBuf(), frq! CASE 6 ' here we use frq! for exponentiation of RND(1) NoiseWave SampBuf(), frq! CASE 7 'gauss GaussWave SampBuf(), frq! CASE 8 'dirac DiracWave SampBuf(), frq! CASE 9 ' expon ExponWave SampBuf(), frq! CASE 10 ' read file-data for sample LeesSampleFile SampBuf(), FileName$ CASE 11 'beta END SELECT SELECT CASE transformkeuze% CASE 1 Samp2Spek SampBuf(), Spectrum() CASE 2 Samp2PowSpek SampBuf(), PwSp() CASE 3 Sm2Pw SampBuf(), PwSp() CASE 4 Sm2PosPw SampBuf(), PwSp() CASE 5 Sm2NegPw SampBuf(), PwSp() CASE 6 Sm2Midi SampBuf(), MuSp() CASE 7 ' no FFT - just waveform synthesis END SELECT SELECT CASE displaykeuze% CASE 1 ShowSpectrum Spectrum() :' complete math with neg. f CASE 2 ShowWave SampBuf() CASE 3 ShowIO SampBuf(), PwSp() CASE 4 ShowPwSp PwSp() CASE 5 ShowReal Spectrum() CASE 6 ShowImag Spectrum() CASE 7 ShowRealImag Spectrum() : ' complete math with negative f's CASE 8 ShowMidiSpectrum SampBuf(), MuSp() CASE 9 ' speed calculation of FFT's CASE 10 ' writes the sample-waveform to one or more files filnam$= "WAV" + RTRIM$(LTRIM$(STR$(frq!))) + ".DTA" SchrijfSampleFile SampBuf(), filnam$ END SELECT count% = count% +1 NEXT frq! IF displaykeuze% = 9 THEN eindetime! = TIMER tijdsduur! = eindetime! - starttime! FFTduur! = (tijdsduur! * 1000) / count% LOCATE 29,10 PRINT "Duration for a single FFT ="; FFTduur!; "ms"; sleep END IF SLEEP END '*************************************************************************** SUB Samp2Spek (SampBuf() AS Single, Spectrum() AS Kompleks) STATIC Pi2!, Toggle% IF Toggle% = 0 THEN Toggle% = -1 Pi2! = 6.283185307 END IF ' deze procedure vormt het hart van de ' Fourier-transformatie. g% = 1 ' = 1 for direct transform ' = -1 for inverse transform NrSamp%= ABS(LBOUND(SampBuf)) + UBOUND(SampBuf) + 1 :' number of samples Dsize% = NrSamp% * 2 ' for each input sample, we need a complex pair ' to calculate the FFT. ' Hence the Size of the internal D-array must be ' twice the amount of samples. ' Of course, if the input data (the samples) is ' complex, the sizes of this array and that of the ' sample-array would be the same. ' This would be the case if we wanted to prove ' the fact that f(t) <=> FFT(f): the Fourier transform ' of a spectrum, is a waveform. DIM D(0 TO Dsize% - 1) AS Single 'D(0),D(2),D(4), D(6)...,D(1022) = reals ' D(1), D(3), D(5), D(7)...,D(1023) = imaginary parts 'D(2046) = SampBuf(-1) 'D(2044) = SampBuf(-2) '... 'D(1024) = SampBuf(LBOUND(SampBuf)) = most negative time 'D(2047) = imaginary part 'D(1025) = most negative time sample imaginary part ' now we copy the sample buffer to the D() array ' in the even places (reals only, conform to ADC cards ' giving real data only) ' therefore D() should have twice as many places ' as SampBuf() ' We put real world samples only in the first half of D() on even places ' not extending beyond the first half of D() ' copy real positive time-samples first: FOR i% = 0 TO UBOUND(SampBuf) D(i%*2) = SampBuf(i%) NEXT i% ' now the real negative time-samples, ' again only on even places: ' this is only used for the analytical cases ' of the transform. IF LBOUND(SampBuf) < 0 THEN FOR i% = -1 TO LBOUND(SampBuf) STEP -1 ' note that i% * 2 is always negative! D(Dsize% + (i% * 2)) = SampBuf(i%) NEXT i% END IF ' We now have D() in a folded cylinder-array ' we go on working now only with D(), since ' the results of the transform are to be found ' in D() at the end of the algorithm ' Step 1:------------------------------------------------------------------ ' Here we unfold the cylinder such that ' the original order: ' 1024 1026 1028... 2047 0 2 4 6 8 10 ... 1022 1023 becomes: ' 0 2 4 ... 1023 1024 1026 1028 1030 1032 1034 ... 2044 2047 ' --negative time date-- positive time data ------------------------- j% = 1 FOR i% = 1 TO Dsize% STEP 2 IF (i% - j%) < 0 THEN SWAP D(j%-1), D(i%-1) :'real part swap: D(0) becomes D(1024) SWAP D(j%), D(i%) :'imag part swap: D(1) becomes D(1025) etc... END IF m% = NrSamp% DO IF (j%-m%) <= 0 THEN EXIT DO j% = j% - m% m% = m% / 2 :'1024-512-256-128-64-32-16-8-4-2-1 IF (m% - 2) < 0 THEN EXIT DO LOOP j% = j% + m% NEXT i% ' Step 2:-----------------Do the Fourier transform------------------------- x% = 2 DO 'IF (x% - Dsize%) >=0 THEN EXIT DO- 'exit condition inserted at the end of the do-loop ' 14.09.97 f% = x% + x% :' 4,8,16,32,64,126,256... h = - Pi2! / (g% * x%) :' g% is 1 or -1 only r = SIN(h/2) w = -2 * r * r v = SIN(h) p = 1 q = 0 FOR m% = 1 TO x% STEP 2 FOR i% = m% TO Dsize% STEP f% :' i% is allways odd j% = i% + x% :' j% is odd t = (p * D(j%-1)) - (q * D(j%)) s = (p * D(j%)) + (q * D(j%-1)) D(j%-1) = D(i%-1) - t :' real D(j%) = D(i%) - s :' imag D(i%-1) = D(i%-1) + t :' real D(i%) = D(i%) + s :' imag NEXT i% t = p p = p + (p * w) - (q * v) q = q + (q * w) + (t * v) NEXT m x% = f% LOOP UNTIL x% >= Dsize% ' Step3: - reorganisation of resulting data--------------------------------- ' refold the array to a closed cilinder... ' warning: 'D(0)-(D((NrSamp%/2)-1)) = positive frequency values ' D(0) - (D(511)) for a 1024 point FFT 'D(1024)- D(1024+511) = negative frequency values!!! ' thus the following would be correct: REDIM Spectrum(-NrSamp%/2 TO (NrSamp%/2)-1) AS Kompleks FOR i% = 0 TO (NrSamp%\2) -1 STEP 2 :' first half of D() - positive frequencies Spectrum(i%\2).real= D(i%) Spectrum(i%\2).imag= D(i%+1) NEXT i% ' now the negative frequencies: pnt% = -1 FOR i% = UBOUND(D) TO (UBOUND(D)/2) STEP -2 :' 2047,2045,2043... 1023 Spectrum(pnt%).real = D(i%-1) Spectrum(pnt%).imag = D(i%) pnt% = pnt% -1 IF pnt% < LBOUND(Spectrum) THEN EXIT FOR NEXT i% END SUB '*************************************************************************** SUB Samp2PowSpek (SampBuf() AS Single, PwSp() AS Single) ' returns a bipolar frequency-spectrum power array ' including negative frequencies. ' For normal periodic transforms the negative half ' will be symmetric with the positive half. REDIM Spectrum(-1 TO 1) AS Kompleks Samp2Spek SampBuf(), Spectrum() REDIM PwSp(LBOUND(Spectrum) TO UBOUND(Spectrum)) AS Single FOR i%= LBOUND(Spectrum) TO UBOUND(Spectrum) PwSp(i%) = SQR((Spectrum(i%).real^2) + (Spectrum(i%).imag^2)) ' for proper scaling we have to do: ' PwSp(i%) = PwSp(i%) /(UBOUND(Spectrum)*2) NEXT i% END SUB '************************************************************************** SUB Sm2Pw (SampBuf() AS Single, PwSp() AS Single) ' experiment REDIM Spectrum(-1 TO 1) AS Kompleks Samp2Spek SampBuf(), Spectrum() REDIM PwSp(0 TO UBOUND(Spectrum)) AS Single ' we sommeren de positieve en negatieve spektraalkomponenten FOR i%= 0 TO UBOUND(Spectrum) frpos =SQR(((Spectrum(i%).real)^2) + ((Spectrum(i%).imag)^2)) frneg =SQR(((Spectrum(-i%).real)^2) + ((Spectrum(-i%).imag)^2)) PwSp(i%) = frpos + frneg ' if you want to rescale do this: 'PwSp(i%) = (frpos + frneg) / (UBOUND(PwSp)*2) NEXT i% END SUB '************************************************************************** SUB Sm2PosPw (SampBuf() AS Single, PwSp() AS Single) REDIM Spectrum(-1 TO 1) AS Kompleks Samp2Spek SampBuf(), Spectrum() REDIM PwSp(0 TO UBOUND(Spectrum)) AS Single 'we return only the positive frequency components... FOR i%= 0 TO UBOUND(Spectrum) frpos =SQR(((Spectrum(i%).real)^2) + ((Spectrum(i%).imag)^2)) PwSp(i%) = frpos NEXT i% END SUB '************************************************************************** SUB Sm2NegPw (SampBuf() AS Single, PwSp() AS Single) REDIM Spectrum(-1 TO 1) AS Kompleks Samp2Spek SampBuf(), Spectrum() REDIM PwSp(0 TO UBOUND(Spectrum)) AS Single 'we return only the negative frequency components... FOR i%= 0 TO UBOUND(Spectrum) frneg = SQR(((Spectrum(-i%).real)^2) + ((Spectrum(-i%).imag)^2)) PwSp(i%) = frneg NEXT i% END SUB '************************************************************************** SUB Sm2Midi (SampBuf() AS Single, MuSp() AS Integer) ' this makes a good chromatic transform... ' the result is usefull as midi-volumes or velocities ' Output seems to be shifted a half tone however... ' We would better shift the tuning of La! a quarter tone ' lower to that the band around the note is returned and not ' the band from a semitone lower to the note itself. STATIC Toggle%, Gronddo! IF Toggle% = 0 THEN La! = 431.000 :' should be quarter tone below La Gronddo! = (La! * (2! ^ (3!/12!))) / 64! Toggle% = -1 END IF ' returns the power spectrum as a midi-array NrSamp% = ABS(LBOUND(SampBuf)) + UBOUND(SampBuf) + 1 REDIM PwSp(0 TO 1) AS Single Sm2PosPw SampBuf(), PwSp() REDIM MuSp(0 TO 127) AS Integer maxval! = 0 oldnoot% = 0 FOR f = 9 TO UBOUND(PwSp) ' no valid midi notes below 9Hz!!! this was a serious bug... noot% = INT(12 * (LOG(f) - LOG(Gronddo!)) / (LOG(2))) ' calculate the sigma sum for the spectral ' power density in a log scale: IF noot% = oldnoot% THEN scv! = scv! + (PwSp(f)/SQR(f)) ELSE scv! = PwSp(f) / SQR(f) oldnoot% = noot% END IF MuSp(noot%)= INT(scv!) IF MuSp(noot%) > maxval! THEN maxval! = MuSp(noot%) NEXT f IF maxval! > 0 THEN coef! = 127 / maxval! FOR i% = 0 TO 127 scv! = MuSp(i%) * coef! MuSp(i%) = INT(scv!) IF MuSp(i%) > 127 THEN MuSp(i%)= 127 IF MuSp(i%) < 0 THEN MuSp(i%)=0 NEXT i% END IF END SUB '******** Waveform generation procedures ********************************** SUB SineWave (SampBuf() AS Single, Freq AS Single) period! = (UBOUND(SampBuf)+1) / Freq :' aantal samples per periode i% = LBOUND(SampBuf) hoekkonstante! = 6.283185307 /period! DO SampBuf(i%) = SIN(hoekkonstante! * i%) INCR i% LOOP UNTIL i% > UBOUND(SampBuf) ' normalized -1 to +1 END SUB '*************************************************************************** SUB SquareWave (SampBuf() AS Single, Freq AS Single) period! = (UBOUND(SampBuf)+ 1) / Freq stap! = 2! / period! staplp! = 1! / stap! i%= LBOUND(SampBuf) DO IF i% <= staplp! THEN SampBuf(i%) = 1 ELSE SampBuf(i%) = -1 IF i% > (staplp! + (1/stap!)) THEN staplp! = staplp! + (2/stap!) END IF END IF INCR i% LOOP UNTIL i% > UBOUND(SampBuf) END SUB '************************************************************************** SUB SawTooth (SampBuf() AS Single, Freq AS Single) period = (UBOUND(SampBuf) + 1) / Freq :' aantal samples per periode stap = 2! / period SampBuf(0)= 0 i% = 1 DO SampBuf(i%)= SampBuf(i%-1) + stap IF SampBuf(i%) > 1 THEN SampBuf(i%)=-1 INCR i% LOOP UNTIL i% > UBOUND(SampBuf) END SUB '************************************************************************** SUB TopHat (SampBuf() AS Single, breedte AS Integer) FOR i% = - breedte/2 TO breedte/2 SampBuf(i%) = 1! / breedte NEXT i% END SUB '************************************************************************** SUB SincWave (SampBuf() AS Single, Freq AS Single) period! = (2 * UBOUND(SampBuf)+1) / Freq :' aantal samples per periode i% = LBOUND(SampBuf) hoekkonstante! = 6.283185307 / period! DO hk! = hoekkonstante! * i% IF hk! = 0 THEN SampBuf(i%) = 1 ELSE SampBuf(i%) = (SIN(hk!)) / hk! END IF INCR i% LOOP UNTIL i% > UBOUND(SampBuf) END SUB '************************************************************************** SUB NoiseWave (SampBuf() AS Single, Freq AS Single) FOR i% = LBOUND(SampBuf) TO UBOUND(SampBuf) tmp! = RND(1) IF tmp! >= 0.5 THEN smp!= -(RND(1) ^ Freq) ELSE smp! = RND(1) ^ Freq END IF SampBuf(i%)= smp! NEXT i% END SUB '************************************************************************** SUB GaussWave (SampBuf() AS Single, Breedte AS Single) ' Freq (here Breedte) is use for the witdh of the Gaussian ' since it counts in units, we have to rescale NrSamp% = ABS(LBOUND(SampBuf)) + UBOUND(SampBuf) + 1 Wd! = Breedte / LOG2(NrSamp%) :' scaling to avoid overflows QWd! = Wd! * Wd! FOR i% = LBOUND(SampBuf) TO UBOUND(SampBuf) it! = (i% * 2!) / NrSamp% Qit! = it! * it! macht# = -(Qit! / Qwd!) gausval! = EXP(macht#) ' = e ^( - (it^2)/(breedte^2)) SampBuf(i%) = gausval! NEXT i% ' the fourier of a gauss-curve is another ' gausscurve! END SUB '************************************************************************** SUB ExponWave (SampBuf() AS Single, Slope AS Single) ' Freq (here Slope) is used for the slope of the exponential ' since it counts in units, we have to rescale NrSamp% = ABS(LBOUND(SampBuf)) + UBOUND(SampBuf) + 1 Wd! = Slope / LOG2(NrSamp%) :' scaling to avoid overflows FOR i% = LBOUND(SampBuf) TO UBOUND(SampBuf) it! = (i% * 2!) / NrSamp% kurval! = EXP(-it!/Wd!) j%= i%/ 2 posidx% = j% + (NrSamp%/4) negidx% = -posidx% IF posidx% <= UBOUND(SampBuf) THEN SampBuf(posidx%) = kurval! END IF IF negidx% >= LBOUND(SampBuf) THEN SampBuf(negidx%) = kurval! END IF NEXT i% END SUB '************************************************************************** SUB LeesSampleFile (SampBuf() AS Single, Fil AS String) OPEN Fil FOR INPUT AS #1 'l% = LEN(#1) FOR i% = LBOUND(SampBuf) TO UBOUND(SampBuf) INPUT #1, SampBuf(i%) NEXT i% ' data is not normalized. If we want normalisation, we have to divide all samples ' by their resolution. CLOSE #1 END SUB '************************************************************************* SUB SchrijfSampleFile (SampBuf() AS Single, Fil AS String) OPEN Fil FOR OUTPUT AS #1 resolution% = 2 ^ 12 FOR i% = LBOUND(SampBuf) TO UBOUND(SampBuf) PRINT #1, INT(SampBuf(i%) * resolution%); NEXT i% ' we suppose the data in SampBuf is normalized. ' The procedure converts to 16-bit signed integer format ' in 12 bit resolution CLOSE #1 END SUB '************************************************************************** SUB DiracWave (SampBuf() AS Single, Points AS Single) FOR i% = LBOUND(SampBuf) TO UBOUND(SampBuf) IF ABS(i%)=Points THEN SampBuf(i%) = 1000000 :' should be infinite... ELSE SampBuf(i%) = 0 END IF NEXT i% ' the fourier transform of a single spike dirac wave is unity! END SUB '******** DISPLAY PROCEDURES ********************************************** SUB ShowIO (SampBuf() AS Single, PwSp() AS Single) SCREEN 12 CLS 1 Size = ABS(LBOUND(SampBuf)) + UBOUND(SampBuf) + 1 LOCATE 1, 1 PRINT "Waveform for"; Size;" samples"; ' upperwindow is for input waveform (sample-array) ' ************************************************ VIEW (20,30)-(620,230),1,13 ' find largest values for autoscaling maxamp = 0 maxtim = LBOUND(SampBuf)-1 FOR i% = LBOUND(SampBuf) TO UBOUND(SampBuf) IF ABS(SampBuf(i%)) > maxamp THEN maxamp= ABS(SampBuf(i%)) maxtim=i% END IF NEXT i% IF maxamp > 0 THEN WINDOW (LBOUND(SampBuf),-maxamp)-(UBOUND(SampBuf),maxamp) ' now display the sample: FOR i% = LBOUND(SampBuf) TO UBOUND(SampBuf) IF i% = 0 THEN LINE (0,-maxamp)-(0,maxamp),15 END IF LINE (i%,0)-(i%,SampBuf(i%)),10 NEXT i END IF ' lowerwindow is for power spectrum display ' ***************************************** LOCATE 16,1 PRINT "Spectral power distribution:"; VIEW (20,260)-(620,460),1,13 ' find largest values maxamp = 0 maxfrq% = LBOUND(PwSp)-1 FOR i% = LBOUND(PwSp) TO UBOUND(PwSp) IF PwSp(i%) > maxamp THEN maxamp = PwSp(i%) maxfrq% = i% END IF NEXT i% LOCATE 16, 40: PRINT "Fmax =";maxfrq%;" Hz @P=";maxamp; IF maxamp > 0 THEN WINDOW (LBOUND(PwSp),0)-(UBOUND(PwSp),maxamp) FOR i% = LBOUND(PwSp) TO UBOUND(PwSp) IF Fraster% (i%) THEN IF i% <> 0 THEN kleur% = 12 ELSE kleur% = 15 LINE (i%,0)-(i%,maxamp),kleur% END IF LINE (i%,0) - (i%,PwSp(i%)),10 NEXT i% END IF BottomLine END SUB '************************************************************************ SUB ShowWave (SampBuf() AS Single) SCREEN 12: CLS 1 VIEW (10,40)-(620,460),1,13 ' find largest values - for autoscaling maxamp = 0 maxtim =LBOUND(SampBuf)- 1 FOR i% = LBOUND(SampBuf) TO UBOUND(SampBuf) IF ABS(SampBuf(i%)) > maxamp THEN maxamp= ABS(SampBuf(i%)) maxtim =i% END IF NEXT i% IF maxamp > 0 THEN WINDOW (LBOUND(SampBuf),-maxamp)-(UBOUND(SampBuf),maxamp) ' now display the sample: FOR i% = LBOUND(SampBuf) TO UBOUND(SampBuf) IF i% = 0 THEN LINE (0,-maxamp)-(0,maxamp),15 :' zero line END IF LINE (i%,0)-(i%,SampBuf(i%)),10 NEXT i% END IF Size = ABS(LBOUND(SampBuf)) + UBOUND(SampBuf) + 1 COLOR 7 LOCATE 1, 1: PRINT "Waveformdisplay over"; Size; "samples Umax=";maxamp; BottomLine END SUB '************************************************************************* SUB ShowSpectrum (Spectrum() AS Kompleks) ' shows the complete mathematic transform, inclusive negative ' frequencies ' Performs the power calculation inside this procedure SCREEN 12 CLS 1 VIEW (10,50)-(620,450),1,13 ' find largest values for autoscaling maxamp = 0 maxfrq =LBOUND(Spectrum)-1 FOR i% = LBOUND(Spectrum) TO UBOUND(Spectrum) power = SQR(((Spectrum(i%).real)^2) + ((Spectrum(i%).imag)^2)) IF power > maxamp THEN maxamp = power maxfrq = i% END IF NEXT i% IF maxamp > 0 THEN WINDOW (LBOUND(Spectrum),0)-(UBOUND(Spectrum),maxamp) ' now display the spectrum: FOR i% = LBOUND(Spectrum) TO UBOUND(Spectrum) IF Fraster% (i%) THEN IF i% <> 0 THEN kleur% = 12 ELSE kleur% = 15 LINE (i%,0)-(i%,maxamp),kleur% END IF power = SQR((Spectrum(i%).real^2) + (Spectrum(i%).imag^2)) LINE (i%,0)-(i%,power),10 NEXT i% END IF Size = ABS(LBOUND(Spectrum)) + UBOUND(Spectrum) + 1 LOCATE 1, 1: PRINT "Spectrumdisplay over"; Size; "points";Size /2; "frequencies"; LOCATE 2, 1: PRINT "f-max=";maxfrq;"Hz","P-max=";maxamp;" "; BottomLine END SUB '*************************************************************************** SUB ShowPwSp (PwSp() AS Single) ' this shows energy in the frequency bands SCREEN 12 CLS 1 VIEW (10,50)-(620,450),1,13 ' find largest values for autoscaling maxamp = 0 maxfrq = LBOUND(PwSp) -1 FOR i% = LBOUND(PwSp) TO UBOUND(PwSp) IF PwSp(i%) > maxamp THEN maxamp = PwSp(i%) maxfrq = i% END IF NEXT i% IF maxamp > 0 THEN WINDOW (LBOUND(PwSp),0)-(UBOUND(PwSp),maxamp) ' now display the spectrum: FOR i% = LBOUND(PwSp) TO UBOUND(PwSp) IF Fraster% (i%) THEN IF i% <> 0 THEN kleur% = 12 ELSE kleur% = 15 LINE (i%,0)-(i%,maxamp),kleur% END IF LINE (i%,0)-(i%,PwSp(i%)),10 NEXT i% END IF LOCATE 1, 1: PRINT "PowerSpectrumdisplay"; LOCATE 2, 1: PRINT "f-max=";maxfrq;"Hz","P-max=";maxamp;" "; BottomLine END SUB '*************************************************************************** SUB ShowMidiSpectrum (SampBuf() AS Single, MuSp() AS Integer) 'dual display procedure: ' upper window = waveform ' lower window = chromatic transform SCREEN 12 CLS 1 NrSamp% = ABS(LBOUND(SampBuf)) + UBOUND(SampBuf) + 1 LOCATE 1, 1 PRINT "Waveform for"; NrSamp%;" samples"; ' upperwindow is for input waveform (sample-array) ' ************************************************ VIEW (20,30)-(620,230),1,13 ' find largest values for autoscaling maxamp = 0 maxtim = LBOUND(SampBuf)-1 FOR i% = LBOUND(SampBuf) TO UBOUND(SampBuf) IF ABS(SampBuf(i%)) > maxamp THEN maxamp= ABS(SampBuf(i%)) maxtim=i% END IF NEXT i% IF maxamp > 0 THEN WINDOW (LBOUND(SampBuf),-maxamp)-(UBOUND(SampBuf),maxamp) ' now display the sample: FOR i% = LBOUND(SampBuf) TO UBOUND(SampBuf) IF i% = 0 THEN LINE (0,-maxamp)-(0,maxamp),15 END IF LINE (i%,0)-(i%,SampBuf(i%)),10 NEXT i END IF ' lowerwindow is for midi spectrum display ' ***************************************** LOCATE 16,1 PRINT "Chromatic Power distribution:"; VIEW (20,260)-(620,460),1,13 ' find largest values maxamp = 0 maxfrq% = LBOUND(MuSp) FOR i% = LBOUND(MuSp) TO UBOUND(MuSp) IF MuSp(i%) > maxamp THEN maxamp = MuSp(i%) maxfrq% = i% END IF NEXT i% LOCATE 16, 40: PRINT "Note-max =";maxfrq%;" Vel=";maxamp; IF maxamp > 0 THEN WINDOW (LBOUND(MuSp),0)-(UBOUND(MuSp),maxamp) FOR i% = LBOUND(MuSp) TO UBOUND(MuSp) IF i% MOD 12 = 0 THEN LINE (i%,0)-(i%,127),4 :' oktaaflijnen END IF SELECT CASE (i% MOD 12) CASE 0,2,4,5,7,9,11 kleur% =15 :' wit CASE ELSE kleur% =7 :' grijs END SELECT LINE (i%,0) - (i%,MuSp(i%)),kleur% NEXT i% END IF BottomLine END SUB '*************************************************************************** SUB ShowReal (Spectrum() AS Kompleks) SCREEN 12 CLS 1 VIEW (10,50)-(620,450),1,13 ' find largest values for autoscaling maxamp = 0 maxfrq = LBOUND(Spectrum) -1 FOR i% = LBOUND(Spectrum) TO UBOUND(Spectrum) p = ABS(Spectrum(i%).real) IF p > maxamp THEN maxamp = p maxfrq = i END IF NEXT i% IF maxamp > 0 THEN WINDOW (LBOUND(Spectrum),-maxamp)-(UBOUND(Spectrum),maxamp) ' now display the real part components of the spectrum: FOR i% = LBOUND(Spectrum) TO UBOUND(Spectrum) IF Fraster% (i%) THEN IF i% <> 0 THEN kleur% = 12 ELSE kleur% = 15 LINE (i%,-maxamp)-(i%,maxamp),kleur% END IF p = Spectrum(i%).real LINE (i%,0)-(i%,p),10 NEXT i% END IF LOCATE 1, 1: PRINT "Real part spectrum display over"; UBOUND(Spectrum)+1; "points"; LOCATE 2, 1: PRINT "f-max=";maxfrq;"Hz"," real-max=";maxamp;" "; BottomLine END SUB '**************************************************************************** SUB ShowImag (Spectrum() AS Kompleks) SCREEN 12 CLS 1 ' autoscaling VIEW (10,50)-(620,450),1,13 ' find largest values maxamp = 0 maxfrq = LBOUND(Spectrum) -1 FOR i% = LBOUND(Spectrum) TO UBOUND(Spectrum) p = ABS(Spectrum(i%).imag) IF p > maxamp THEN maxamp = p maxfrq = i% END IF NEXT i% IF maxamp > 0 THEN WINDOW (LBOUND(Spectrum),-maxamp)-(UBOUND(Spectrum),maxamp) ' now display the spectrum: FOR i% = LBOUND(Spectrum) TO UBOUND(Spectrum) IF Fraster% (i%) THEN IF i% <> 0 THEN kleur% = 12 ELSE kleur% = 15 LINE (i%,-maxamp)-(i%,maxamp),kleur% END IF p = Spectrum(i%).imag LINE (i%,0)-(i%,p),10 NEXT i% END IF LOCATE 1, 1: PRINT "Imaginary part spectrum display over"; UBOUND(Spectrum)+1; "points"; LOCATE 2, 1: PRINT "f-max=";maxfrq;"Hz"," imag-max=";maxamp;" "; BottomLine END SUB '**************************************************************************** SUB ShowRealImag (Spectrum() AS Kompleks) ' also shows negative frequency components! SCREEN 12 CLS 1 ' upper window = real number part - autoscaling VIEW (10,50)-(620,200),1,13 ' find largest values maxreal = 0 maximag = 0 maxfr = LBOUND(Spectrum)-1 maxfi = LBOUND(Spectrum)-1 ' find both extremes , but we set the scale identical for both ' windows to allow easy comparisons. FOR i% = LBOUND(Spectrum) TO UBOUND(Spectrum) p = ABS(Spectrum(i%).real) pv = ABS(Spectrum(i%).imag) IF p > maxreal THEN maxreal = p maxfr = i% END IF IF pv > maximag THEN maximag = pv maxfi = i% END IF NEXT i% IF maxreal > maximag THEN maxamp=maxreal ELSE maxamp= maximag IF maxamp > 0 THEN WINDOW (LBOUND(Spectrum),-maxamp)-(UBOUND(Spectrum),maxamp) ' now display the real part components of the spectrum: FOR i% = LBOUND(Spectrum) TO UBOUND(Spectrum) IF Fraster% (i%) THEN IF i% <> 0 THEN kleur% = 12 ELSE kleur% = 15 LINE (i%,-maxamp)-(i%,maxamp),kleur% END IF p = Spectrum(i%).real LINE (i%,0)-(i%,p),10 NEXT i END IF LOCATE 1, 1: PRINT "Real part spectrum display over"; (UBOUND(Spectrum)+1)*2; "points ";UBOUND(Spectrum);" frequencies"; LOCATE 2, 1: PRINT "f-max=";maxfr;"Hz","real-max=";maxreal;" "; ' lower window = imaginary numbers VIEW (10,250)-(620,400),1,13 IF maxamp > 0 THEN WINDOW (LBOUND(Spectrum),-maxamp)-(UBOUND(Spectrum),maxamp) FOR i% = LBOUND(Spectrum) TO UBOUND(Spectrum) IF Fraster% (i%) THEN IF i% <> 0 THEN kleur% = 12 ELSE kleur% = 15 LINE (i%,-maxamp)-(i%,maxamp),kleur% END IF p = Spectrum(i%).imag LINE (i%,0)-(i%,p),10 NEXT i% END IF LOCATE 28, 1: PRINT "Imaginary part spectrum display over"; (UBOUND(Spectrum)+1)*2; "points"; LOCATE 29, 1: PRINT "f-max=";maxfi;"Hz","imag-max=";maximag;" "; BottomLine END SUB '*************************************************************************** ' utility procedures FUNCTION Fraster% (i AS Integer) ' calculates the frequency points on the x-axis for the graphs on powers of 2 Fraster% = 0 IF i = 0 THEN Fraster% = -1 EXIT FUNCTION END IF FOR j% = 2 TO 14 IF ABS(i) = 2 ^ j% THEN Fraster% = -1 : EXIT FOR NEXT j% END FUNCTION SUB BottomLine COLOR 1 LOCATE 30,15: PRINT "FFT-explorer by Prof.Dr.Godfried-Willem Raes"; COLOR 7 END SUB ' E.O.F.