mrandri19/fourier.html

## fourier.html
<!DOCTYPE html>
<html>
<head>
	<meta charset="utf-8">
	<title>DFT</title>
	<style>
		h1 {
			font-family: sans-serif;
		}
		#canvas_container {
			display: flex;
			flex-wrap: wrap;
			justify-content: flex-start;
		}

		#canvas_container > * {
			margin: 0.2rem;
		}
	</style>
</head>
<body>
	<h1>Discrete Fourier Transform</h1>
	<div id="canvas_container"></div>

	<script src="fourier.js"></script>
</body>
</html>

## fourier.js
//*****************************************************************************
function make_function(fn, N) {
	const x = new Float32Array(N)

	for(let n = 0; n < N; n++) {
		x[n] = fn(n)
	}

	return x
}

function make_sin(f, N) {
	return make_function(n => Math.sin(2 * Math.PI * f * n / N), N)
}

function make_exp(a, N) {
	return make_function(n => Math.exp(- a * n / N), N)
}

function make_sinc(f, N) {
	return make_function(
		n => (n == 0) ? 1 : Math.sin(
			2 * Math.PI * f * n / N) / (2 * Math.PI * f * n / N),
			N
		)
}

function make_delta(N) {
	return make_function(n => (n == 0) ? 1 : 0, N)
}

function make_gaussian(mu, sigma, N) {
	return make_function(
		n => Math.exp(
			- 0.5 * Math.pow((n - mu) / sigma, 2)
			) / (sigma * Math.sqrt(2 * Math.PI)),
		N
		)
}

function make_rect(l, N) {
	return make_function(n => (n < l) ? 1 : 0, N)
}


//*****************************************************************************
function delay(x, d) {
	const N = x.length

	const y = new Float32Array(N)

	for(let n = 0; n < d; n++) {
		y[n] = 0
	}
	for(let n = d; n < N; n++) {
		y[n] = x[n - d]
	}

	return y
}

function invert(x) {
	const N = x.length

	const y = new Float32Array(N)

	for(let n = 0; n < N; n++) {
		y[n] = x[N - n - 1]
	}

	return y
}

function add(...xs) {
	const N = xs[0].length

	if (!xs.map(x => x.length == N).reduce((a, b) => a && b)) {
		throw "All signals should have the same length"
	}

	const z = new Float32Array(N)
	for(let n = 0; n < N; n++) {
			z[n] = xs.map(x => x[n]).reduce((a, b) => a + b)
	}
	return z
}

//*****************************************************************************
function magnitude(X_re, X_im) {
	const N = X_re.length

	const X_mag = new Float32Array(N)

	for(let k = 0; k < N; k++) {
		X_mag[k] = Math.sqrt(X_re[k] * X_re[k] + X_im[k] * X_im[k])
	}

	return X_mag
}

function phase(X_re, X_im) {
	const N = X_re.length

	const X_phase = new Float32Array(N)

	for(let k = 0; k < N; k++) {
		X_phase[k] = Math.atan(X_im[k], X_re[k])
	}

	return X_phase
}

//*****************************************************************************
function cross_correlation_finite_power(x, y) {
	const N = x.length

	const R = new Float32Array(N)

	for(let n = 0; n < N; n++) {
		for(let k = 0; k < (N - n); k++) {
			R[n] += (x[k + n] * y[k])
		}

		R[n] /= N
	}

	return R
}

//*****************************************************************************
function DFT(x) {
	const N = x.length

	const X_re = new Float32Array(N)
	const X_im = new Float32Array(N)

	for(let k = 0; k < N; k++) {
		for(let n = 0; n < N; n++) {
			X_re[k] += (x[n] * Math.cos(- 2 * Math.PI * n * k / N))
			X_im[k] += (x[n] * Math.sin(- 2 * Math.PI * n * k / N))
		}
	}

	return [X_re, X_im]
}

function FFT_recursive(x) {
	const N = x.length

	if(Math.log2(N) % 1 !== 0) {
		throw "The signal's length must be a power of 2"
	}

	const e = new Float32Array(N/2)
	const o = new Float32Array(N/2)

	for(let k = 0; k < N/2; k++) {
		e[k] = x[2 * k]
		o[k] = x[2 * k + 1]
	}

	const [E_re, E_im] = (N == 2) ? DFT(e) : FFT_recursive(e)
	const [O_re, O_im] = (N == 2) ? DFT(o) : FFT_recursive(o)

	const X_re = new Float32Array(N)
	const X_im = new Float32Array(N)

	for(let k = 0; k < N/2; k++) {
		const t_re = Math.cos(- 2 * Math.PI * 1 * k / N)
		const t_im = Math.sin(- 2 * Math.PI * 1 * k / N)

		const O_t_re = O_re[k] * t_re - O_im[k] * t_im
		const O_t_im = O_im[k] * t_re + O_re[k] * t_im

		X_re[k] = E_re[k] + O_t_re
		X_im[k] = E_im[k] + O_t_im

		X_re[k + N/2] = E_re[k] - O_t_re
		X_im[k + N/2] = E_im[k] - O_t_im
	}

	return [X_re, X_im]
}

//*****************************************************************************
function fftshift(x) {
	const N = x.length
	if(Math.log2(N) % 1 !== 0) {
		throw "The signal's length must be a power of 2"
	}

	const y = new Float32Array(N)

	for(let n = 0; n < N; n++) {
		y[(n + (N / 2)) % N] = x[n]
	}

	return y
}

//*****************************************************************************
function main() {
	const N =  Math.pow(2, 14)
	console.log(N)

	{
		const x = add(
			make_sin(10, N),
			make_sin(100, N),
			make_sin(200, N),
			// delay(make_sinc(50, N), N/2),
			// invert(delay(make_sinc(50, N), N-N/2)),
			// make_exp(50, N),
			// delay(make_delta(N), 512),
			// invert(delay(make_gaussian(0, 10, N), 512)),
			// delay(make_gaussian(0, 10, N), 512),
			// delay(make_rect(50,N),100),
		)

		{
			// TODO: handle case where N != K, i.e. where I want the K-point DFT
			// on a N-point signal

			console.time('FFT_recursive')
			const [X_re, X_im] = FFT_recursive(x)
			console.timeEnd('FFT_recursive')

			const X_mag = magnitude(X_re, X_im)
			plot(fftshift(X_mag), "|X[k]| (FFT)")
		}

		// {
		// 	console.time('DFT')
		// 	const [X_re, X_im] = DFT(x)
		// 	console.timeEnd('DFT')

		// 	const X_mag = magnitude(X_re, X_im)
		// 	plot(fftshift(X_mag), "|X[k]| (DFT)")
		// }

	}
}

function plot(x, title) {
	const N = x.length

	const canvas_container = document.getElementById("canvas_container")
	const canvas = document.createElement("canvas")
	canvas_container.appendChild(canvas)

	const W = canvas.width = 500
	const H = canvas.height = 3/4 * W

	const c = canvas.getContext("2d")

	// draw background
	c.fillStyle = "#eee"
	c.fillRect(0, 0, W, H)


	// normalization constant
	const max_x = x.reduce((a, b) => (a > b) ? a : b)

	// draw signal
	c.fillStyle = "#000"

	const baseline = H / 2
	const max_amplitude = H / 3
	for(let n = 0; n < N; n++) {
		x_norm = x[n] / max_x
		c.fillRect(
			n / N * W,
			 baseline - (x_norm * max_amplitude),
			1,
			(x_norm * max_amplitude)
		)
	}

	// x axis
	c.strokeStyle = "#00f"
	c.beginPath()
	c.moveTo(0, H / 2)
	c.lineTo(W, H / 2)
	c.stroke()

	c.font = "24px sans-serif"
	c.textAlign = "center"
	c.fillText(title, W / 2, 24)

}

//*****************************************************************************
main()
	<!DOCTYPE html>
	<html>
	<head>
	<meta charset="utf-8">
	<title>DFT</title>
	<style>
	h1 {
	font-family: sans-serif;
	}
	#canvas_container {
	display: flex;
	flex-wrap: wrap;
	justify-content: flex-start;
	}

	#canvas_container > * {
	margin: 0.2rem;
	}
	</style>
	</head>
	<body>
	<h1>Discrete Fourier Transform</h1>
	<div id="canvas_container"></div>

	<script src="fourier.js"></script>
	</body>
	</html>
	//*****************************************************************************
	function make_function(fn, N) {
	const x = new Float32Array(N)

	for(let n = 0; n < N; n++) {
	x[n] = fn(n)
	}

	return x
	}

	function make_sin(f, N) {
	return make_function(n => Math.sin(2 * Math.PI * f * n / N), N)
	}

	function make_exp(a, N) {
	return make_function(n => Math.exp(- a * n / N), N)
	}

	function make_sinc(f, N) {
	return make_function(
	n => (n == 0) ? 1 : Math.sin(
	2 * Math.PI * f * n / N) / (2 * Math.PI * f * n / N),
	N
	)
	}

	function make_delta(N) {
	return make_function(n => (n == 0) ? 1 : 0, N)
	}

	function make_gaussian(mu, sigma, N) {
	return make_function(
	n => Math.exp(
	- 0.5 * Math.pow((n - mu) / sigma, 2)
	) / (sigma * Math.sqrt(2 * Math.PI)),
	N
	)
	}

	function make_rect(l, N) {
	return make_function(n => (n < l) ? 1 : 0, N)
	}


	//*****************************************************************************
	function delay(x, d) {
	const N = x.length

	const y = new Float32Array(N)

	for(let n = 0; n < d; n++) {
	y[n] = 0
	}
	for(let n = d; n < N; n++) {
	y[n] = x[n - d]
	}

	return y
	}

	function invert(x) {
	const N = x.length

	const y = new Float32Array(N)

	for(let n = 0; n < N; n++) {
	y[n] = x[N - n - 1]
	}

	return y
	}

	function add(...xs) {
	const N = xs[0].length

	if (!xs.map(x => x.length == N).reduce((a, b) => a && b)) {
	throw "All signals should have the same length"
	}

	const z = new Float32Array(N)
	for(let n = 0; n < N; n++) {
	z[n] = xs.map(x => x[n]).reduce((a, b) => a + b)
	}
	return z
	}

	//*****************************************************************************
	function magnitude(X_re, X_im) {
	const N = X_re.length

	const X_mag = new Float32Array(N)

	for(let k = 0; k < N; k++) {
	X_mag[k] = Math.sqrt(X_re[k] * X_re[k] + X_im[k] * X_im[k])
	}

	return X_mag
	}

	function phase(X_re, X_im) {
	const N = X_re.length

	const X_phase = new Float32Array(N)

	for(let k = 0; k < N; k++) {
	X_phase[k] = Math.atan(X_im[k], X_re[k])
	}

	return X_phase
	}

	//*****************************************************************************
	function cross_correlation_finite_power(x, y) {
	const N = x.length

	const R = new Float32Array(N)

	for(let n = 0; n < N; n++) {
	for(let k = 0; k < (N - n); k++) {
	R[n] += (x[k + n] * y[k])
	}

	R[n] /= N
	}

	return R
	}

	//*****************************************************************************
	function DFT(x) {
	const N = x.length

	const X_re = new Float32Array(N)
	const X_im = new Float32Array(N)

	for(let k = 0; k < N; k++) {
	for(let n = 0; n < N; n++) {
	X_re[k] += (x[n] * Math.cos(- 2 * Math.PI * n * k / N))
	X_im[k] += (x[n] * Math.sin(- 2 * Math.PI * n * k / N))
	}
	}

	return [X_re, X_im]
	}

	function FFT_recursive(x) {
	const N = x.length

	if(Math.log2(N) % 1 !== 0) {
	throw "The signal's length must be a power of 2"
	}

	const e = new Float32Array(N/2)
	const o = new Float32Array(N/2)

	for(let k = 0; k < N/2; k++) {
	e[k] = x[2 * k]
	o[k] = x[2 * k + 1]
	}

	const [E_re, E_im] = (N == 2) ? DFT(e) : FFT_recursive(e)
	const [O_re, O_im] = (N == 2) ? DFT(o) : FFT_recursive(o)

	const X_re = new Float32Array(N)
	const X_im = new Float32Array(N)

	for(let k = 0; k < N/2; k++) {
	const t_re = Math.cos(- 2 * Math.PI * 1 * k / N)
	const t_im = Math.sin(- 2 * Math.PI * 1 * k / N)

	const O_t_re = O_re[k] * t_re - O_im[k] * t_im
	const O_t_im = O_im[k] * t_re + O_re[k] * t_im

	X_re[k] = E_re[k] + O_t_re
	X_im[k] = E_im[k] + O_t_im

	X_re[k + N/2] = E_re[k] - O_t_re
	X_im[k + N/2] = E_im[k] - O_t_im
	}

	return [X_re, X_im]
	}

	//*****************************************************************************
	function fftshift(x) {
	const N = x.length
	if(Math.log2(N) % 1 !== 0) {
	throw "The signal's length must be a power of 2"
	}

	const y = new Float32Array(N)

	for(let n = 0; n < N; n++) {
	y[(n + (N / 2)) % N] = x[n]
	}

	return y
	}

	//*****************************************************************************
	function main() {
	const N = Math.pow(2, 14)
	console.log(N)

	{
	const x = add(
	make_sin(10, N),
	make_sin(100, N),
	make_sin(200, N),
	// delay(make_sinc(50, N), N/2),
	// invert(delay(make_sinc(50, N), N-N/2)),
	// make_exp(50, N),
	// delay(make_delta(N), 512),
	// invert(delay(make_gaussian(0, 10, N), 512)),
	// delay(make_gaussian(0, 10, N), 512),
	// delay(make_rect(50,N),100),
	)

	{
	// TODO: handle case where N != K, i.e. where I want the K-point DFT
	// on a N-point signal

	console.time('FFT_recursive')
	const [X_re, X_im] = FFT_recursive(x)
	console.timeEnd('FFT_recursive')

	const X_mag = magnitude(X_re, X_im)
	plot(fftshift(X_mag), "\|X[k]\| (FFT)")
	}

	// {
	// console.time('DFT')
	// const [X_re, X_im] = DFT(x)
	// console.timeEnd('DFT')

	// const X_mag = magnitude(X_re, X_im)
	// plot(fftshift(X_mag), "\|X[k]\| (DFT)")
	// }

	}
	}

	function plot(x, title) {
	const N = x.length

	const canvas_container = document.getElementById("canvas_container")
	const canvas = document.createElement("canvas")
	canvas_container.appendChild(canvas)

	const W = canvas.width = 500
	const H = canvas.height = 3/4 * W

	const c = canvas.getContext("2d")

	// draw background
	c.fillStyle = "#eee"
	c.fillRect(0, 0, W, H)


	// normalization constant
	const max_x = x.reduce((a, b) => (a > b) ? a : b)

	// draw signal
	c.fillStyle = "#000"

	const baseline = H / 2
	const max_amplitude = H / 3
	for(let n = 0; n < N; n++) {
	x_norm = x[n] / max_x
	c.fillRect(
	n / N * W,
	baseline - (x_norm * max_amplitude),
	1,
	(x_norm * max_amplitude)
	)
	}

	// x axis
	c.strokeStyle = "#00f"
	c.beginPath()
	c.moveTo(0, H / 2)
	c.lineTo(W, H / 2)
	c.stroke()

	c.font = "24px sans-serif"
	c.textAlign = "center"
	c.fillText(title, W / 2, 24)

	}

	//*****************************************************************************
	main()