Discrete Sine Transform via the FFT

Mon, Feb 5, 2018

Several months ago, I had a quick look at a recent paper describing how to use Wavelets to price options under stochastic volatility models with a known characteristic function. The more classic method is to use some numerical quadrature directly on the Fourier integral as described in this paper for example. When I read the paper, I was skeptical about the Wavelet approach, since it looked complicated, and with many additional parameters.

I recently took a closer look and it turns out my judgment was bad. The additional parameters can be easily automatically and reliably set by very simple rules that work well (a subsequent paper by the same author, which applies the method to Bermudan options, clarifies this). The method is also not as complex as I first imagined. And more importantly, the FFT makes it fast. It is quite amazing to see the power of the FFT in action. It really is because of the FFT that the Shannon Wavelet method is practical.

Anyway one of the things that need to be computed are the payoff coefficients, and one expression is just the sum of a discrete Cosine transform (DCT) and a discrete Sine transform (DST). I was wondering then about a simple way to use the FFT for the Sine transform. There are many papers around how to use the FFT to compute the Cosine transform. A technique that is efficient and simple is the one of Makhoul.

The coefficients that need to be computed for all k can be represented by the following equation $$ V_{k} = \sum_{j=0}^{N-1} a_j \cos\left(\pi k\frac{j+\frac{1}{2}}{N}\right) + b_j \sin\left(\pi k\frac{j+\frac{1}{2}}{N}\right) $$ with $ N= 2^{\bar{J}-1} $ for some positive integer $ \bar{J} $. Makhoul algorithm to compute the DCT of size N with one FFT of size N consists in

initialize the FFT coefficients $c_j$ with: $$ c_j = a_{2j} \quad,,\quad c_{N-1-j} = a_{2j+1} \quad \text{ for } j = 0,…, \frac{N}{2}-1 $$
and then from the result of the FFT $ hat{c} $, the DCT coefficients $ \hat{a} $ are $$ \hat{a}_k = \Re\left[ \hat{c}_j e^{-i \pi\frac{k}{2N}} \right],. $$

Makhoul does not specify the equivalent formula for the DST, but we can do something similar.

We first initialize the FFT coefficients $ c_j $ with: $$ c_j = b_{2j} \quad,,\quad c_{N-1-j} = -b_{2j+1} \quad \text{ for } j = 0,…, \frac{N}{2}-1 $$
and then from the result of the FFT $\hat{c}$, the DST coefficients $\hat{b}$ are $$\hat{b}_k = -\Im\left[ \hat{c}_j e^{-i \pi\frac{k}{2N}} \right],. $$

For maximum performance, the two FFTs can reuse the same sine and cosine tables. And the last step of the DCT and DST can be combined together.

Another approach would be to compute a single FFT of size 2N as we can rewrite the coefficients as $$ V_{k} = \Re\left[e^{-i \pi \frac{k}{2N}} \sum_{j=0}^{2N-1} (a_j + i b_j) e^{-2i\pi k\frac{j}{2N}} \right] $$ with $a_j = b_j =0 $ for $ j \geq N $

In fact the two are almost equivalent, since a FFT of size 2N is decomposed internally into two FFTs of size N.

With this in mind we can improve a bit the above to merge the two transforms together:

We first initialize the FFT coefficients $ c_j $ with: $$ c_j = a_{2j} + i b_{2j} \quad,,\quad c_{N-1-j} = a_{2j+1} -i b_{2j+1} \quad \text{ for } j = 0,…, \frac{N}{2}-1 $$
and then from the result of the FFT $\hat{c}$, the coefficients are $$V_k = \Re\left[ \hat{c}_j e^{-i \pi\frac{k}{2N}} \right],. $$

And we have computed the coefficients with a single FFT of size N.

In terms of Go language code, the classic type-2 sine Fourier transforms reads

func DST2(vector []float64) error {
	n := len(vector)
	result := make([]float64, n)
	for k := range vector {
		sum := 0.0
		for j, vj := range vector {
			sum += math.Sin((float64(j)+0.5)*math.Pi*float64(k)/float64(n)) * vj
		}
		result[k] = sum
	}
	copy(vector, result)
	return nil
}

And the Fast Type-2 Sine transform reads

func FastDST2ByFFT(vector []float64) error {
	n := len(vector)
	halfLen := n / 2
	z := make([]complex128, n)
	for i := 0; i < halfLen; i++ {
		z[i] = complex(vector[i*2], 0)
		z[n-1-i] = complex(-vector[i*2+1], 0)
	}
	if n&1 == 1 {
		z[halfLen] = complex(vector[n-1], 0)
	}
	err := TransformRadix2Complex(z)
	for i := range z {
		temp := float64(i) * math.Pi / float64(n*2)
		s, c := math.Sincos(temp)
		vector[i] = real(z[i])*s - imag(z[i])*c
	}
	return err
}

The function TransformRadix2Complex, used above, steams from a standard implementation of the discrete Fourier transform (DFT) of the given complex vector, storing the result back into the vector (The vector’s length must be a power of 2. Uses the Cooley-Tukey decimation-in-time radix-2 algorithm).

func TransformRadix2Complex(z []complex128) error {
	n := len(z)
	levels := uint(31 - bits.LeadingZeros32(uint32(n))) // Equal to floor(log2(n))
	if (1 << levels) != n {
		return fmt.Errorf("Length is not a power of 2")
	}

	// Trigonometric tables
	cosTable, sinTable := MakeCosSinTable(n)

	// Bit-reversed addressing permutation
	for i := range z {
		j := int(bits.Reverse32(uint32(i)) >> (32 - levels))
		if j > i {
			z[i], z[j] = z[j], z[i]
		}
	}

	// Cooley-Tukey decimation-in-time radix-2 FFT
	for size := 2; size <= n; size *= 2 {
		halfsize := size / 2
		tablestep := n / size
		for i := 0; i < n; i += size {
			j := i
			k := 0
			for ; j < i+halfsize; j++ {
				l := j + halfsize
				tpre := real(z[l])*cosTable[k] + imag(z[l])*sinTable[k]
				tpim := -real(z[l])*sinTable[k] + imag(z[l])*cosTable[k]
				z[l] = z[j] - complex(tpre, tpim)
				z[j] += complex(tpre, tpim)
				k += tablestep
			}
		}
		if size == n { // Prevent overflow in 'size *= 2'
			break
		}
	}

	return nil
}

func MakeCosSinTable(n int) ([]float64, []float64) {
	cosTable := make([]float64, n/2)
	sinTable := make([]float64, n/2)
	for i := range cosTable {
		s, c := math.Sincos(2 * math.Pi * float64(i) / float64(n))
		cosTable[i] = c
		sinTable[i] = s
	}
	return cosTable, sinTable
}

Discrete Sine Transform via the FFT

Comments