diff --git a/noise/perlin.go b/noise/perlin.go
new file mode 100644
index 0000000..07cd30f
--- /dev/null
+++ b/noise/perlin.go
@@ -0,0 +1,151 @@
+package noise
+
+import (
+	"math"
+	"math/rand"
+)
+
+const perlinB = 0x100
+const perlinBM = 0xff
+const perlinN = 0x1000
+
+// Perlin contains a perlin noise function (both 1D and 2D).
+type Perlin struct {
+	Seed int64
+	p    [perlinB + perlinB + 2]int
+	g1   [perlinB + perlinB + 2]float64
+	g2   [perlinB + perlinB + 2][2]float64
+}
+
+// NewPerlin creates a new perlin noise function with seed seed.
+func NewPerlin(seed int64) *Perlin {
+	p := &Perlin{Seed: seed}
+	p.Init()
+	return p
+}
+
+func perlinSmoothCurve(t float64) float64 { return t * t * (3. - 2.*t) }
+
+func perlinLerp(t, a, b float64) float64 { return a + t*(b-a) }
+
+func perlinAt2D(q [2]float64, rx, ry float64) float64 { return rx*q[0] + ry*q[1] }
+
+func (p *Perlin) noise1d(x float64) float64 {
+	tx := x + perlinN
+	bx0 := int(tx) & perlinBM
+	bx1 := (bx0 + 1) & perlinBM
+	rx0 := tx - math.Floor(tx)
+	rx1 := rx0 - 1.
+
+	sx := perlinSmoothCurve(rx0)
+	u := rx0 * p.g1[p.p[bx0]]
+	v := rx1 * p.g1[p.p[bx1]]
+
+	return perlinLerp(sx, u, v)
+}
+
+func (p *Perlin) noise2d(x, y float64) float64 {
+	tx := x + perlinN
+	bx0 := int(tx) & perlinBM
+	bx1 := (bx0 + 1) & perlinBM
+	rx0 := tx - math.Floor(tx)
+	rx1 := rx0 - 1.
+
+	ty := y + perlinN
+	by0 := int(ty) & perlinBM
+	by1 := (by0 + 1) & perlinBM
+	ry0 := ty - math.Floor(ty)
+	ry1 := ry0 - 1.
+
+	i := p.p[bx0]
+	j := p.p[bx1]
+
+	b00 := p.p[i+by0]
+	b10 := p.p[j+by0]
+	b01 := p.p[i+by1]
+	b11 := p.p[j+by1]
+
+	sx := perlinSmoothCurve(rx0)
+	sy := perlinSmoothCurve(ry0)
+
+	u := perlinAt2D(p.g2[b00], rx0, ry0)
+	v := perlinAt2D(p.g2[b10], rx1, ry0)
+	a := perlinLerp(sx, u, v)
+
+	u = perlinAt2D(p.g2[b01], rx0, ry1)
+	v = perlinAt2D(p.g2[b11], rx1, ry1)
+	b := perlinLerp(sx, u, v)
+
+	return perlinLerp(sy, a, b)
+}
+
+func normalize2D(v [2]float64) [2]float64 {
+	var s float64
+
+	s = math.Sqrt(v[0]*v[0] + v[1]*v[1])
+	v[0] /= s
+	v[1] /= s
+	return v
+}
+
+// Init initialised the Perlin noise source. Must be called unless struct was created with NewPerlin(...).
+func (p *Perlin) Init() {
+	r := rand.New(rand.NewSource(p.Seed))
+
+	for i := 0; i < perlinB; i++ {
+		p.p[i] = i
+		p.g1[i] = float64((r.Int()%(perlinB+perlinB))-perlinB) / perlinB
+
+		for j := 0; j < 2; j++ {
+			p.g2[i][j] = float64((r.Int()%(perlinB+perlinB))-perlinB) / perlinB
+		}
+		p.g2[i] = normalize2D(p.g2[i])
+		for j := 0; j < 3; j++ {
+			r.Int()
+		}
+	}
+
+	for i := perlinB; i > 0; i-- {
+		j := r.Int() % perlinB
+		p.p[i], p.p[j] = p.p[j], p.p[i]
+	}
+
+	for i := 0; i < perlinB+2; i++ {
+		p.p[perlinB+i] = p.p[i]
+		p.g1[perlinB+i] = p.g1[i]
+		for j := 0; j < 2; j++ {
+			p.g2[perlinB+i][j] = p.g2[i][j]
+		}
+	}
+}
+
+// Noise1D generates the noise value for x. alpha defines the weight when the sum is formed. Typically it is 2 (function is noisier when alpha approaches 1). beta  is the harmonic scaling/spacing (typically 2). n defines the number of harmonics.
+func (p *Perlin) Noise1D(x, alpha, beta float64, n int) float64 {
+	var val, sum float64
+	var scale float64 = 1
+	var coord float64 = x
+
+	for i := 0; i < n; i++ {
+		val = p.noise1d(coord)
+		sum += val / scale
+		scale *= alpha
+		coord *= beta
+	}
+	return sum
+}
+
+// Noise2D generates the noise value for coordinate (x, y). alpha defines the weight when the sum is formed. Typically it is 2 (function is noisier when alpha approaches 1). beta  is the harmonic scaling/spacing (typically 2). n defines the number of harmonics.
+func (p *Perlin) Noise2D(x, y, alpha, beta float64, n int) float64 {
+	var val, sum float64
+	var scale float64 = 1
+	var coord = [2]float64{x, y}
+
+	for i := 0; i < n; i++ {
+		val = p.noise2d(coord[0], coord[1])
+		sum += val / scale
+		scale *= alpha
+		coord[0] *= beta
+		coord[1] *= beta
+	}
+	return sum
+}