Module:Harmonic entropy: Difference between revisions

← Older edit

@@ Line 1: / Line 1: @@
-local getArgs = require('Module:Arguments').getArgs
+local limits = require("Module:Limits")
-local u = require('Module:Utils')
+local rat = require("Module:Rational")
 local p = {}
--- return measure in cents of an interval ratio, rounded to prec decimal places
+-- Compute Harmonic Shannon entropy for an interval of `c` cents
-function p.to_cents(frame)
+-- `c`, `deviation`: in cents
-	local args = getArgs(frame)
+-- `ratios`: an array of rational numbers
-	return p._to_cents(args[1], args[2])
+-- `norm`: a function of rational numbers
-end
+function p.harmonic_entropy(c, ratios, deviation, norm)
+	norm = norm or function(ratio)
+		return math.sqrt(rat.benedetti_height(ratio))
+	end
+	deviation = deviation or 1200 * math.log(1.01, 2)
+	ratios = ratios
+		or limits.integer_limit(200, function(ratio)
+			if math.abs(rat.cents(ratio) - c) > 3 * deviation then
+				return 1 / 0
+			end
+			return norm(ratio)
+		end, 100)
+	local function gaussian(x)
+		return math.exp(-x * x / (2 * deviation * deviation)) / (deviation * math.sqrt(2 * math.pi))
+	end
+	local function weighted_gaussian(ratio)
+		return gaussian(rat.cents(ratio) - c) / norm(ratio)
+	end
+	local q_norm = 0
+	for _, ratio in pairs(ratios) do
+		q_norm = q_norm + weighted_gaussian(ratio)
+	end
+	local function probability(ratio)
+		return weighted_gaussian(ratio) / q_norm
+	end
-function p._to_cents(ratio, prec)
+	local entropy = 0
-	-- ratio defaults to 1
+	for _, ratio in pairs(ratios) do
-	ratio = u.eval_num_arg(ratio, 1)
+		local p_i = probability(ratio)
-	-- output without rounding if prec is omitted
+		if p_i > 1e-5 then
-	if prec then
+			entropy = entropy - p_i * math.log(p_i)
-		return u._round(1200*u._log(ratio), prec)
+		end
-	-- output with rounding if prec is given
-	else
-		-- prec defaults to 3
-		prec = u.eval_num_arg(prec, 3)
-		return 1200*u._log(ratio)
 	end
+	return entropy
 end
 return p