PROGRAM analyze
! determine the Fourier coefficients a_k and b_k
CALL parameters(N,nmax,delta,period)
CALL screen(nmax,period,#1,#2)
CALL coefficents(N,nmax,delta,period,#1,#2)
END

SUB parameters(N,nmax,delta,period)
    INPUT prompt "number of data points N (even) = ": N
    INPUT prompt "sampling time dt = ": delta
    LET period = N*delta          ! assumed period
    ! maximum value of mode corresponding to Nyquist frequency
    LET nmax = N/2
END SUB

SUB screen(nmax,period,#1,#2)
    LET ymax = 2
    LET ticksize = ymax/50
    OPEN #1: screen 0,1,0.5,1
    PRINT "     a_k";
    PRINT "        ";
    PRINT using "frequency interval = #.#####": 2*pi/period
    SET WINDOW -1,nmax+1,-ymax,ymax
    CALL plotaxis(nmax,ticksize)
    SET COLOR "red"
    OPEN #2: screen 0,1,0,0.5
    PRINT "     b_k"
    SET WINDOW -1,nmax+1,-ymax,ymax
    CALL plotaxis(nmax,ticksize)
    SET COLOR "red"
END SUB

SUB plotaxis(nmax,ticksize)
    PLOT LINES: 0,0;nmax,0
    FOR k = 1 to nmax
        PLOT LINES: k,-ticksize;k,ticksize
    NEXT k
END SUB

SUB coefficents(N,nmax,delta,period,#1,#2)
    DECLARE DEF f
    FOR k = 0 to nmax
        LET ak = 0
        LET bk = 0
        LET wk = 2*pi*k/period
        ! rectangular approximation
        FOR i = 0 to N - 1
            LET t = i*delta
            LET ak = ak + f(t)*cos(wk*t)
            LET bk = bk + f(t)*sin(wk*t)
        NEXT i
        LET ak = 2*ak/N
        LET bk = 2*bk/N
        WINDOW #1
        PLOT LINES: k,0;k,ak
        WINDOW #2
        PLOT LINES: k,0;k,bk
    NEXT k
END SUB

FUNCTION f(t)
    LET w0 = 0.1*pi
    LET f = sin(w0*t)             ! simple example
END DEF