Skip to content

Instantly share code, notes, and snippets.

@mrandri19

mrandri19/fourier.html

Created January 23, 2020 22:10
Show Gist options
  • Save mrandri19/d0f2716ea221ea277377f0a987a63a92 to your computer and use it in GitHub Desktop.
Save mrandri19/d0f2716ea221ea277377f0a987a63a92 to your computer and use it in GitHub Desktop.
Download ZIP
Experiments with the Discrete Fourier Transform
Raw
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>
Raw
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()
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment