' frequency modulator & demodulator for fm_radio_template.reason
' by veezay
'
' FModem.exe (fmdemod9.exe)
' latest version 2018-11-27 (prev 2018-07-31)
'
' new in version:
' - doesn't overwrite output file fmodem_out.wav, instead uses running numbering from 001 to 999
'
'
' stuff todo:
' - add timestamp of start & finishing ETA of processing
' - select new RANDOMIZE seed from input audio
' - rewrite WAV header as a TYPE
'
'
' stuff to try out:
' - figure out a nice carrier frequency? or dunno? maybe the current is just fine? depends on stuff.
' - try different carrier waveforms such as triangle (cleanest zero crossing calculation?)


RANDOMIZE TIMER


DIM lutsize AS INTEGER, lutsizeoffset AS INTEGER, lfolutsize AS INTEGER, audiolutsize AS INTEGER, lutindex AS _UNSIGNED INTEGER, lutindextemp AS DOUBLE
DIM palutratio AS _UNSIGNED LONG, lfopalutratio AS _UNSIGNED LONG, audiopalutratio AS _UNSIGNED LONG
DIM outputsample(2) AS INTEGER, previousoutputsample(2) AS INTEGER
DIM sinelookup(1024, 1024) AS INTEGER
DIM modulatoroutput AS INTEGER
DIM samplerate AS _UNSIGNED LONG, newsample AS INTEGER, samplex AS _UNSIGNED LONG
DIM carrier AS _UNSIGNED LONG, carrieros AS _UNSIGNED LONG
DIM currentsample AS _UNSIGNED _INTEGER64, s AS _UNSIGNED _INTEGER64, position AS _UNSIGNED _INTEGER64, lastposition AS _UNSIGNED _INTEGER64
DIM sampleoffset AS _UNSIGNED _INTEGER64, outputoffset AS _UNSIGNED _INTEGER64, samplebuffer AS _UNSIGNED LONG, outputbuffer AS _UNSIGNED LONG, n AS _UNSIGNED LONG
DIM dummy AS _BYTE, dummyint AS INTEGER
DIM downsampled AS LONG, writedownsampled AS _BIT, lastchunk AS _BIT
lastchunk = 0

oversampling = 4 ' 8
samplerate = 88200 * oversampling
carrier = 45059 '35059 '* 2 + 1
carrieros = carrier \ oversampling
writedownsampled = 1

' lookup table sizes
lfolutsize = 64 ' not in use currently
audiolutsize = 256
lfopalutratio = 2 ^ 32 / lfolutsize ' phase accumulator size / lookup table size ratio

' how many input samples to load per chunk + same for output
samplebuffer = 2 ^ 16
outputbuffer = 2 ^ 16
outputfileextra = 256 ' compensation for me being a shitty coder

' make lookup tables of oscillator
GOSUB generatelookuptable


' define input file
IF COMMAND$ = "" THEN
    inputfiles$ = "C:\files\reason\exports\fm_in.wav" ' probably should replace this with some error message telling user to plz drag&drop .wav on .exe
ELSE
    inputfiles$ = COMMAND$
    inputfile = 1
    FOR i = 1 TO LEN(inputfiles$)
        char$ = MID$(inputfiles$, i, 1)
        IF char$ = " " THEN
            inputfile = inputfile + 1
        ELSE
            inputfile$(inputfile) = inputfile$(inputfile) + char$
        END IF
    NEXT i
    numberofinputs = inputfile
    FOR i = 1 TO numberofinputs
        LOCATE i: PRINT " input file "; trim$(i); ": "; inputfile$(i)
    NEXT i
END IF

' start timing the execution
DIM time AS DOUBLE
time = TIMER(0.001)

' get input file specifics
DIM maxchannels AS INTEGER, numberofchannels(numberofinputs) AS INTEGER, samplelength(numberofinputs) AS _UNSIGNED _INTEGER64, sampletime(numberofinputs)
DIM inputfilenumber(numberofinputs)

FOR i = 1 TO numberofinputs
    inputfilenumber(i) = FREEFILE
    OPEN inputfile$(i) FOR BINARY AS inputfilenumber(i)
    GET inputfilenumber(i), 23, numberofchannels(i)

    IF maxchannels < numberofchannels(i) THEN maxchannels = numberofchannels(i)

    samplelength(i) = (LOF(inputfilenumber(i)) - 44) / 2 / numberofchannels(i)
    sampletime(i) = samplelength(i) / 88200 ' we assume our input sample rate is always 88200 Hz, which it of course always is
NEXT i

' initialize output file
DIM zcrossprev(maxchannels) AS _FLOAT, zcrosscurr(maxchannels) AS _FLOAT, samplecount(maxchannels) AS _UNSIGNED LONG, writtensamples(maxchannels) AS _UNSIGNED LONG
DIM sample(numberofinputs, samplebuffer * oversampling, maxchannels) AS INTEGER
IF writedownsampled THEN
    DIM outputfile(outputbuffer * maxchannels / oversampling) AS INTEGER ' the actual buffer to be written to file
ELSE
    DIM outputfile(outputbuffer * maxchannels) AS INTEGER ' the actual buffer to be written to file
END IF
DIM outputfilecopy((outputbuffer + outputfileextra) * maxchannels) AS INTEGER ' the buffer that is being filled + shifted
DIM newinputsample(i, maxchannels) AS INTEGER, oldinputsample(i, maxchannels) AS INTEGER
DIM bufferfull(maxchannels) AS INTEGER

' write output to the same directory as the first input file
FOR i = LEN(inputfile$(1)) TO 1 STEP -1
    IF MID$(inputfile$(1), i, 1) = "\" THEN EXIT FOR
NEXT i
outputfile$ = MID$(inputfile$(1), 1, i) + "fmodem_out_001.wav"

filenameid = 1

' let's see if fmodem_out_001.wav exists and increment the number in case that happens to be the case
WHILE _FILEEXISTS(outputfile$)
    'KILL outputfile$

    filenameoffset = LEN(outputfile$) - 6
    outputfile$ = MID$(outputfile$, 0, filenameoffset)
    filenameid = filenameid + 1
    IF filenameid = 1000 THEN
        PRINT " too many outputs, i'm broken. remove fmodem_out_xxx.wavs" ' :D
        END
    END IF

    filenameidstr$ = ""
    IF filenameid < 10 THEN filenameidstr$ = "0"
    IF filenameid < 100 THEN filenameidstr$ = filenameidstr$ + "0"
    filenameidstr$ = filenameidstr$ + LTRIM$(RTRIM$(STR$(filenameid)))
    outputfile$ = outputfile$ + filenameidstr$ + ".wav"
WEND

outputfilenumber = FREEFILE
OPEN outputfile$ FOR BINARY AS outputfilenumber

' rewrite this shit into TYPE
DIM foolong AS LONG, fooint AS INTEGER
foostr$ = "RIFF": PUT outputfilenumber, 1, foostr$
foostr$ = "WAVE": PUT outputfilenumber, 9, foostr$
foostr$ = "fmt ": PUT outputfilenumber, 13, foostr$
foolong = 16: PUT outputfilenumber, 17, foolong
fooint = 1: PUT outputfilenumber, 21, fooint
PUT outputfilenumber, 23, maxchannels
foolong = samplerate
IF writedownsampled THEN foolong = foolong / oversampling
PUT outputfilenumber, 25, foolong
foolong = samplerate * 2 * maxchannels
IF writedownsampled THEN foolong = foolong / oversampling
PUT outputfilenumber, 29, foolong
fooint = 2 * maxchannels: PUT outputfilenumber, 33, fooint
fooint = 16: PUT outputfilenumber, 35, fooint
foostr$ = "data": PUT outputfilenumber, 37, foostr$


LOCATE 2 + numberofinputs: PRINT " making progress"
progresstimer = _FREETIMER
ON TIMER(progresstimer, 0.1) GOSUB printprogress
TIMER(progresstimer) ON


' initialize carrier gain specifics
DIM carriergain(numberofinputs) AS DOUBLE, initial_carriergain(numberofinputs) AS DOUBLE, carrierlimitmin(numberofinputs) AS DOUBLE, carrierlimitmax(numberofinputs) AS DOUBLE
' range 0 -> 1 (logarithmic), 0 being no signal at all, 1 good strong signal
' magic numbers defining main signal strength
carrierlimitmin(1) = 0.00035
carrierlimitmax(1) = 0.01120 '0.01280
carriergain(1) = RND * (carrierlimitmax(1) - carrierlimitmin(1)) + carrierlimitmin(1)
' give rest of input signals lower gain limits
FOR i = 2 TO numberofinputs
    carrierlimitmin(i) = 0.00001
    carrierlimitmax(i) = 0.00060
    carriergain(i) = RND * (carrierlimitmax(i) - carrierlimitmin(i)) + carrierlimitmin(i)
NEXT i

' update the carrier signal gains drift every 2^nth sample
DIM carrierdrift_counter AS _UNSIGNED _BIT * 6
'initial_carriergain = carriergain


DIM phaseacc(numberofinputs, maxchannels) AS _UNSIGNED LONG, frequency(numberofinputs) AS _UNSIGNED LONG

FOR i = 1 TO numberofinputs
    frequency(i) = carrier / samplerate * 2 ^ 32
NEXT i

DIM gain AS LONG, phasenoisegain(numberofinputs) AS LONG, phasenoisegainhalf(numberofinputs) AS LONG, noisegain AS SINGLE, noisegainhalf AS SINGLE, demodulationgain AS DOUBLE
gain = 8192 ' so this this is obviously the bandwidth of the carrier signal)
phasenoisegain(1) = 150 ^ 3 ' 200 ^ 3 seems to give a fair bit of noise
phasenoisegainhalf(1) = phasenoisegain(1) / 2
noisegain = 10: noisegainhalf = noisegain / 2
demodulationgain = 32768 / gain ' plz

' lfos not used, also in case need lfos: rewrite using arrays
lfo1 = 2
lfo2 = 3

lfo1rate = 10000
'lfo1gain = 1
lfo2rate = 10000
'lfo1gain = 1

'phaseacc(lfo1) = 0


sampleoffset = 0

' initialize zero crossings
FOR channel = 0 TO maxchannels - 1
    samplecount(channel) = 0
    writtensamples(channel) = 0
    zcrosscurr(channel) = 0 ' initial zero crossing "happened" at the first sample
NEXT channel


' we'll just figure these out here in case we don't need to recalculate them in the loop, which we don't
lutsize = audiolutsize '+ modulatoroutput ' lfo modulating audio lookup table size (should probably interpolate if ever use it)
audiopalutratio = 4294967296 / lutsize ' phase accumulator size / lookup table size ratio

outputmodulatedsignal = 0
IF outputmodulatedsignal = 1 THEN
    IF _FILEEXISTS("signal.fm") THEN KILL "signal.fm"
    OPEN "signal.fm" FOR BINARY AS 10
END IF

carrierdrift_counter = 1 ' start from 1 in case we don't want to mess about with the signal gain

DIM random_drift(numberofinputs) AS DOUBLE, random_drift_multiplier AS DOUBLE
random_drift_multiplier = 0.0005
' the following calculation keeps the "speed" of the occasional interference noise somewhat constant
' regardless of oversampling setting or number of input signals used by scaling the multiplier
' according to how many RND calls are being made
random_drift_multiplier = random_drift_multiplier / oversampling / (2 * numberofinputs + 1)


' main loop
FOR s = 0 TO samplelength(1) * oversampling - 1
    IF s MOD samplebuffer = 0 THEN GOSUB loadsamplechunk

    carrierdrift_counter = carrierdrift_counter + 1 ' comment this out for testing without drift
    IF carrierdrift_counter = 0 THEN
        FOR i = 1 TO numberofinputs
            random_drift(i) = (RND - 0.5) * random_drift_multiplier
            carriergain(i) = carriergain(i) + random_drift(i) ' linear gain drift seems to work fine (also it's less expensive than logarithmic so win/win)
            'carriergain = carriergain + (initial_carriergain - carriergain) * 0.000015 ' autotune (no, not that kind)
            IF carriergain(i) < carrierlimitmin(i) THEN carriergain(i) = carrierlimitmin(i) ' hard clip the gain values
            IF carriergain(i) > carrierlimitmax(i) THEN carriergain(i) = carrierlimitmax(i) ' -||-
        NEXT i
    END IF

    FOR channel = 0 TO maxchannels - 1
        samplecount(channel) = samplecount(channel) + 1

        currentsample = s - sampleoffset + samplebuffer

        outputsample(channel) = 0
        FOR i = 1 TO numberofinputs
            ' following two lines are not necessary to compute in the loop since NOTHING'S MODULATING
            'lutsize = audiolutsize '+ modulatoroutput ' lfo modulating audio lookup table size (should probably interpolate if ever use it)
            'audiopalutratio = 4294967296 / lutsize ' phase accumulator size / lookup table size ratio
            phaseacc(i, channel) = phaseacc(i, channel) + frequency(i) + sample(i, currentsample, channel) * gain + phasenoisegain(i) * RND - phasenoisegainhalf(i)
            audiolutindex(i) = phaseacc(i, channel) / audiopalutratio

            outputsample(channel) = outputsample(channel) + sinelookup(lutsize, audiolutindex(i)) * carriergain(i) '+ noisegain(i) * RND - noisegainhalf(i)
        NEXT i

        ' add noise to the fm'd signal
        outputsample(channel) = outputsample(channel) + noisegain * RND - noisegainhalf

        IF SGN(previousoutputsample(channel)) <> SGN(outputsample(channel)) THEN GOSUB countzerocrossing

        previousoutputsample(channel) = outputsample(channel)

    NEXT channel

    ' debug, output fm'd signal
    IF outputmodulatedsignal THEN PUT #10, , outputsample(0)

    IF INP(96) = 1 THEN
        prematureexit = 1
        EXIT FOR
    END IF
NEXT s


IF prematureexit = 0 THEN

    lastchunk = 1
    GOSUB writeoutputchunk

    foolong = LOF(outputfilenumber) - 8: PUT outputfilenumber, 5, foolong
    foolong = LOF(outputfilenumber) - 44: PUT outputfilenumber, 41, foolong
    CLOSE ' both input and output files
    CLS
    TIMER(progresstimer) OFF

    FOR i = 1 TO numberofinputs
        LOCATE i: PRINT " input file "; trim$(i); ": "; inputfile$(i)
    NEXT i
    LOCATE 2 + numberofinputs: PRINT " finished everything"
    LOCATE 3 + numberofinputs: PRINT " 100.00 % ("; outputfile$; ")"

    IF TIMER(.001) < time THEN time = time - 86400
    executiontime = INT((TIMER(.001) - time) * 100) / 100
    executiontoinputratio = INT(executiontime / sampletime(1) * 1000) / 1000
    LOCATE 5 + numberofinputs: PRINT " execution time:"; executiontime; "seconds (";
    IF executiontoinputratio < 1 THEN PRINT "0";
    PRINT RTRIM$(LTRIM$(STR$(executiontoinputratio))); "x)"

    FOR i = 1 TO numberofinputs
        LOCATE 7 + numberofinputs + i: PRINT " gain "; trim$(i); ": ";
        PRINT USING "#.#########"; carriergain(i)
    NEXT i
ELSE
    CLOSE ' both/all input and output files
    'IF _FILEEXISTS(outputfile$) THEN KILL outputfile$
    CLS
    TIMER(progresstimer) OFF
    LOCATE 1: PRINT " exited."
END IF


END


countzerocrossing:
' calculate the instantaneous frequency from two adjacent zero crossings with linear interpolation
zcrossprev(channel) = zcrosscurr(channel)
zcrosscurr(channel) = s + (previousoutputsample(channel) / (previousoutputsample(channel) - outputsample(channel)))
newsample = (samplerate / (2 * oversampling) / (zcrosscurr(channel) - zcrossprev(channel)) - carrieros) * demodulationgain

FOR n = 1 TO samplecount(channel)
    position = writtensamples(channel) * maxchannels + channel - outputoffset
    outputfilecopy(position) = newsample
    writtensamples(channel) = writtensamples(channel) + 1
NEXT n

IF position >= outputbuffer THEN bufferfull(channel) = 1
bcount = 0
FOR n = 0 TO maxchannels - 1
    bcount = bcount + bufferfull(n)
NEXT n
IF bcount = maxchannels THEN
    GOSUB writeoutputchunk
    'PRINT outputoffset;
END IF

samplecount(channel) = 0
RETURN


printprogress:
progress$ = LTRIM$(RTRIM$(STR$(FIX(writtensamples(0) / (samplelength(1) * oversampling) * 10000))))
WHILE LEN(progress$) <= 2
    progress$ = "0" + progress$
WEND
progress$ = STRING$(5 - LEN(progress$), " ") + MID$(progress$, 1, LEN(progress$) - 2) + "." + MID$(progress$, LEN(progress$) - 1, 2) + " %"
LOCATE 3 + numberofinputs: PRINT progress$
RETURN



loadsamplechunk:
FOR i = 1 TO numberofinputs
    GET inputfilenumber(i), 45 + sampleoffset * numberofchannels(i) * 2 \ oversampling, dummybyte

    FOR samplex = 0 TO samplebuffer \ oversampling - 1
        FOR channel = 0 TO numberofchannels(i) - 1
            GET inputfilenumber(i), , newinputsample(i, channel)
            FOR os = 0 TO oversampling - 1
                ' linear interpolation
                sample(i, samplex * oversampling + os, channel) = oldinputsample(i, channel) - (oldinputsample(i, channel) - newinputsample(i, channel)) * ((os + 1) / oversampling)
            NEXT os
            oldinputsample(i, channel) = newinputsample(i, channel)
        NEXT channel
    NEXT samplex
NEXT i

sampleoffset = sampleoffset + samplebuffer
RETURN


writeoutputchunk:
IF writedownsampled THEN
    ' sample averaging downsampling method
    FOR n = 0 TO outputbuffer - 1 STEP (oversampling * maxchannels)
        FOR c = 0 TO maxchannels - 1
            downsampled = 0
            FOR m = 0 TO oversampling - 1
                ' since sample data is interlaced for all channels, we need to skip in the array to get the samples of a given channel
                downsampled = downsampled + outputfilecopy(n + c + m * maxchannels)
            NEXT m
            downsampled = downsampled / oversampling
            outputfile(n / oversampling + c) = downsampled 'outputfilecopy(n)
        NEXT c
    NEXT n
    PUT outputfilenumber, 45 + (outputoffset * 2 / oversampling), outputfile()
ELSE
    ' output oversampled data as is
    FOR n = 0 TO outputbuffer - 1
        outputfile(n) = outputfilecopy(n)
    NEXT n
    PUT outputfilenumber, 45 + (outputoffset * 2), outputfile()
END IF
IF lastchunk THEN RETURN

FOR n = 0 TO outputfileextra - 1
    outputfilecopy(n) = outputfilecopy(outputbuffer + n)
NEXT n
FOR n = outputfileextra TO outputbuffer
    outputfilecopy(n) = 0
NEXT n

FOR n = 0 TO maxchannels - 1
    bufferfull(n) = 0
NEXT n

outputoffset = outputoffset + outputbuffer
RETURN



generatelookuptable:
' look-up tables containing the sinewave period
FOR size = 1 TO 1023
    FOR i = 0 TO size - 1
        sinelookup(size, i) = -SIN(i / size * _PI * 2) * 16384
    NEXT i
NEXT size
RETURN

FUNCTION trim$ (i)
trim$ = LTRIM$(RTRIM$(STR$(i)))
END FUNCTION