Fourier sine transform from Numerical Recipes

sinft.jl

function four1!(data, nn, isign)
    n = Int(2 * nn)
    j = 1

    for i = 1:2:n
        if j > i
            tempr = data[j]
            tempi = data[j + 1]
            data[j] = data[i]
            data[j + 1] = data[i + 1]
            data[i] = tempr
            data[i + 1] = tempi
        end
        m = n ÷ 2
        while m >= 2 && j > m
            j -= m
            m = m ÷ 2
        end
        j += m
    end
    mmax = 2
    while n > mmax
        istep = 2 * mmax
        theta = 2 * pi / (isign * mmax)
        wpr = -2.0 * sin(0.50 * theta)^2
        wpi = sin(theta)
        wr = 1.0
        wi = 0.0
        for m = 1:2:mmax
            for i = m:istep:n
                j = i + mmax
                tempr = wr * data[j] - wi * data[j + 1]
                tempi = wr * data[j + 1] + wi * data[j]
                data[j] = data[i] - tempr
                data[j + 1] = data[i + 1] - tempi
                data[i] = data[i] + tempr
                data[i + 1] = data[i + 1] + tempi
            end
            wtemp = wr
            wr = wr * wpr - wi * wpi + wr
            wi = wi * wpr + wtemp * wpi + wi
        end
        mmax = istep
    end
end

function realft!(data, n, isign)
    theta = pi / n
    c1 = 0.5
    if isign == 1
        c2 = -0.5
        four1!(data, n, 1)
    else
        c2 = 0.5
        theta = -theta
    end
    wpr = -2.0 * sin(0.5 * theta)^2
    wpi = sin(theta)
    wr = 1.0 + wpr
    wi = wpi
    n2p3 = 2 * n + 3

    for i = 2:(n ÷ 2)
        i1 = 2 * i - 1
        i2 = i1 + 1
        i3 = n2p3 - i2
        i4 = i3 + 1
        wrs = wr
        wis = wi
        h1r = c1 * (data[i1] + data[i3])
        h1i = c1 * (data[i2] - data[i4])
        h2r = -c2 * (data[i2] + data[i4])
        h2i = c2 * (data[i1] - data[i3])
        data[i1] = h1r + wrs * h2r - wis * h2i
        data[i2] = h1i + wrs * h2i + wis * h2r
        data[i3] = h1r - wrs * h2r + wis * h2i
        data[i4] = -h1i + wrs * h2i + wis * h2r
        wtemp = wr
        wr = wr * wpr - wi * wpi + wr
        wi = wi * wpr + wtemp * wpi + wi
    end
    if isign == 1
        h1r = data[1]
        data[1] = h1r + data[2]
        data[2] = h1r - data[2]
    else
        h1r = data[1]
        data[1] = c1 * (h1r + data[2])
        data[2] = c1 * (h1r - data[2])
        four1!(data, n, -1)
    end
end

function sinft!(y, n)
    theta = pi / n
    wr = 1.0
    wi = 0.0
    wpr = -2.0 * sin(0.5 * theta)^2
    wpi = sin(theta)
    y[1] = 0.0
    m = n ÷ 2
    for j = 1:m
        wtemp = wr
        wr = wr * wpr - wi * wpi + wr
        wi = wi * wpr + wtemp * wpi + wi
        y1 = wi * (y[j + 1] + y[n - j + 1])
        y2 = 0.5 * (y[j + 1] - y[n - j + 1])
        y[j + 1] = y1 + y2
        y[n - j + 1] = y1 - y2
    end
    realft!(y, m, 1)
    sum = 0.0
    y[1] *= 0.5
    y[2] = 0.0
    for j = 1:2:(n - 1)
        sum = sum + y[j]
        y[j] = y[j + 1]
        y[j + 1] = sum
    end
end