diff --git a/src/ipa/libipa/awb_bayes.cpp b/src/ipa/libipa/awb_bayes.cpp
new file mode 100644
index 000000000000..8ab6ed661a3e
--- /dev/null
+++ b/src/ipa/libipa/awb_bayes.cpp
@@ -0,0 +1,456 @@
+/* SPDX-License-Identifier: LGPL-2.1-or-later */
+/*
+ * Copyright (C) 2019, Raspberry Pi Ltd
+ * Copyright (C) 2024 Ideas on Board Oy
+ *
+ * Implementation of a bayesian AWB algorithm
+ */
+
+#include "awb_bayes.h"
+
+#include <cmath>
+
+#include <libcamera/base/log.h>
+#include <libcamera/control_ids.h>
+
+#include "colours.h"
+
+/**
+ * \file awb_bayes.h
+ * \brief Implementation of bayesian auto white balance algorithm
+ *
+ * This implementation is based on the initial implementation done by
+ * RaspberryPi.
+ * \todo: Documentation
+ *
+ * \todo Not all the features implemented by RaspberryPi were ported over to
+ * this algorithm because they either rely on hardware features not generally
+ * available or were considered not important enough at the moment.
+ *
+ * The following parts are not implemented:
+ *
+ * - min_pixels: minimum proportion of pixels counted within AWB region for it
+ *   to be "useful"
+ * - min_g: minimum G value of those pixels, to be regarded a "useful"
+ * - min_regions: number of AWB regions that must be "useful" in order to do the
+ *   AWB calculation
+ * - deltaLimit: clamp on colour error term (so as not to penalize non-grey
+ *   excessively)
+ * - bias_proportion: The biasProportion parameter adds a small proportion of
+ *   the counted pixels to a region biased to the biasCT colour temperature.
+ *   A typical value for biasProportion would be between 0.05 to 0.1.
+ * - bias_ct: CT target for the search bias
+ * - sensitivityR: red sensitivity ratio (set to canonical sensor's R/G divided
+ *   by this sensor's R/G)
+ * - sensitivityB: blue sensitivity ratio (set to canonical sensor's B/G divided
+ *   by this sensor's B/G)
+ */
+
+namespace libcamera {
+
+LOG_DECLARE_CATEGORY(Awb)
+
+namespace ipa {
+
+/**
+ * \brief Step size control for CT search
+ */
+constexpr double kSearchStep = 0.2;
+
+/**
+ * \copydoc Interpolator::interpolate()
+ */
+template<>
+void Interpolator<Pwl>::interpolate(const Pwl &a, const Pwl &b, Pwl &dest, double lambda)
+{
+	dest = Pwl::combine(a, b,
+			    [=](double /*x*/, double y0, double y1) -> double {
+				    return y0 * (1.0 - lambda) + y1 * lambda;
+			    });
+}
+
+/**
+ * \class AwbBayes
+ * \brief Implementation of a bayesian auto white balance algorithm
+ *
+ * In a bayesian AWB algorithm the auto white balance estimation is improved by
+ * taking the likelihood of a given lightsource based on the estimated lux level
+ * into account. E.g. If it is very bright we can assume that we are outside and
+ * that colour temperatures around 6500 are preferred.
+ *
+ * The second part of this algorithm is the search for the most likely colour
+ * temperature. It is implemented in AwbBayes::coarseSearch() and in
+ * AwbBayes::fineSearch(). The search works very well without prior likelihoods
+ * and therefore the algorithm itself provides very good results even without
+ * prior likelihoods.
+ */
+
+/**
+ * \var AwbBayes::transversePos_
+ * \brief How far to wander off CT curve towards "more purple"
+ */
+
+/**
+ * \var AwbBayes::transverseNeg_
+ * \brief How far to wander off CT curve towards "more green"
+ */
+
+/**
+ * \var AwbBayes::currentMode_
+ * \brief The currently selected mode
+ */
+
+int AwbBayes::init(const YamlObject &tuningData)
+{
+	int ret = colourGainCurve_.readYaml(tuningData["colourGains"], "ct", "gains");
+	if (ret) {
+		LOG(Awb, Error)
+			<< "Failed to parse 'colourGains' "
+			<< "parameter from tuning file";
+		return ret;
+	}
+
+	ctR_.clear();
+	ctB_.clear();
+	for (const auto &[ct, g] : colourGainCurve_.data()) {
+		ctR_.append(ct, 1.0 / g[0]);
+		ctB_.append(ct, 1.0 / g[1]);
+	}
+
+	/* We will want the inverse functions of these too. */
+	ctRInverse_ = ctR_.inverse().first;
+	ctBInverse_ = ctB_.inverse().first;
+
+	ret = readPriors(tuningData);
+	if (ret) {
+		LOG(Awb, Error) << "Failed to read priors";
+		return ret;
+	}
+
+	ret = parseModeConfigs(tuningData, controls::AwbAuto);
+	if (ret) {
+		LOG(Awb, Error)
+			<< "Failed to parse mode parameter from tuning file";
+		return ret;
+	}
+	currentMode_ = &modes_[controls::AwbAuto];
+
+	transversePos_ = tuningData["transversePos"].get<double>(0.01);
+	transverseNeg_ = tuningData["transverseNeg"].get<double>(0.01);
+	if (transversePos_ <= 0 || transverseNeg_ <= 0) {
+		LOG(Awb, Error) << "AwbConfig: transversePos/Neg must be > 0";
+		return -EINVAL;
+	}
+
+	return 0;
+}
+
+int AwbBayes::readPriors(const YamlObject &tuningData)
+{
+	const auto &priorsList = tuningData["priors"];
+	std::map<uint32_t, Pwl> priors;
+
+	if (!priorsList) {
+		LOG(Awb, Error) << "No priors specified";
+		return -EINVAL;
+	}
+
+	for (const auto &p : priorsList.asList()) {
+		if (!p.contains("lux")) {
+			LOG(Awb, Error) << "Missing lux value";
+			return -EINVAL;
+		}
+
+		uint32_t lux = p["lux"].get<uint32_t>(0);
+		if (priors.count(lux)) {
+			LOG(Awb, Error) << "Duplicate prior for lux value " << lux;
+			return -EINVAL;
+		}
+
+		std::vector<uint32_t> temperatures =
+			p["ct"].getList<uint32_t>().value_or(std::vector<uint32_t>{});
+		std::vector<double> probabilities =
+			p["probability"].getList<double>().value_or(std::vector<double>{});
+
+		if (temperatures.size() != probabilities.size()) {
+			LOG(Awb, Error)
+				<< "Ct and probability array sizes are unequal";
+			return -EINVAL;
+		}
+
+		if (temperatures.empty()) {
+			LOG(Awb, Error)
+				<< "Ct and probability arrays are empty";
+			return -EINVAL;
+		}
+
+		std::map<int, double> ctToProbability;
+		for (unsigned int i = 0; i < temperatures.size(); i++) {
+			int t = temperatures[i];
+			if (ctToProbability.count(t)) {
+				LOG(Awb, Error) << "Duplicate ct value";
+				return -EINVAL;
+			}
+
+			ctToProbability[t] = probabilities[i];
+		}
+
+		auto &pwl = priors[lux];
+		for (const auto &[ct, prob] : ctToProbability) {
+			pwl.append(ct, prob);
+		}
+	}
+
+	if (priors.empty()) {
+		LOG(Awb, Error) << "No priors specified";
+		return -EINVAL;
+	}
+
+	priors_.setData(std::move(priors));
+
+	return 0;
+}
+
+void AwbBayes::handleControls(const ControlList &controls)
+{
+	auto mode = controls.get(controls::AwbMode);
+	if (mode) {
+		auto it = modes_.find(static_cast<controls::AwbModeEnum>(*mode));
+		if (it != modes_.end())
+			currentMode_ = &it->second;
+		else
+			LOG(Awb, Error) << "Unsupported AWB mode " << *mode;
+	}
+}
+
+RGB<double> AwbBayes::gainsFromColourTemperature(double colourTemperature)
+{
+	/*
+	 * \todo: In the RaspberryPi code, the ct curve was interpolated in
+	 * the white point space (1/x) not in gains space. This feels counter
+	 * intuitive, as the gains are in linear space. But I can't prove it.
+	 */
+	const auto &gains = colourGainCurve_.getInterpolated(colourTemperature);
+	return { { gains[0], 1.0, gains[1] } };
+}
+
+AwbResult AwbBayes::calculateAwb(const AwbStats &stats, int lux)
+{
+	ipa::Pwl prior;
+	if (lux > 0) {
+		prior = priors_.getInterpolated(lux);
+		prior.map([](double x, double y) {
+			LOG(Awb, Debug) << "(" << x << "," << y << ")";
+		});
+	} else {
+		prior.append(0, 1.0);
+	}
+
+	double t = coarseSearch(prior, stats);
+	double r = ctR_.eval(t);
+	double b = ctB_.eval(t);
+	LOG(Awb, Debug)
+		<< "After coarse search: r " << r << " b " << b << " (gains r "
+		<< 1 / r << " b " << 1 / b << ")";
+
+	/*
+	 * Original comment from RaspberryPi:
+	 * Not entirely sure how to handle the fine search yet. Mostly the
+	 * estimated CT is already good enough, but the fine search allows us to
+	 * wander transversely off the CT curve. Under some illuminants, where
+	 * there may be more or less green light, this may prove beneficial,
+	 * though I probably need more real datasets before deciding exactly how
+	 * this should be controlled and tuned.
+	 */
+	fineSearch(t, r, b, prior, stats);
+	LOG(Awb, Debug)
+		<< "After fine search: r " << r << " b " << b << " (gains r "
+		<< 1 / r << " b " << 1 / b << ")";
+
+	return { { { 1.0 / r, 1.0, 1.0 / b } }, t };
+}
+
+double AwbBayes::coarseSearch(const ipa::Pwl &prior, const AwbStats &stats) const
+{
+	std::vector<Pwl::Point> points;
+	size_t bestPoint = 0;
+	double t = currentMode_->ctLo;
+	int spanR = -1;
+	int spanB = -1;
+
+	/* Step down the CT curve evaluating log likelihood. */
+	while (true) {
+		double r = ctR_.eval(t, &spanR);
+		double b = ctB_.eval(t, &spanB);
+		RGB<double> gains({ 1 / r, 1.0, 1 / b });
+		double delta2Sum = stats.computeColourError(gains);
+		double priorLogLikelihood = prior.eval(prior.domain().clamp(t));
+		double finalLogLikelihood = delta2Sum - priorLogLikelihood;
+
+		LOG(Awb, Debug) << "Coarse search t: " << t
+				<< " gains: " << gains
+				<< " error: " << delta2Sum
+				<< " prior: " << priorLogLikelihood
+				<< " likelihood: " << finalLogLikelihood;
+
+		points.push_back({ { t, finalLogLikelihood } });
+		if (points.back().y() < points[bestPoint].y())
+			bestPoint = points.size() - 1;
+
+		if (t == currentMode_->ctHi)
+			break;
+
+		/*
+		 * Ensure even steps along the r/b curve by scaling them by the
+		 * current t.
+		 */
+		t = std::min(t + t / 10 * kSearchStep, currentMode_->ctHi);
+	}
+
+	t = points[bestPoint].x();
+	LOG(Awb, Debug) << "Coarse search found CT " << t;
+
+	/*
+	 * We have the best point of the search, but refine it with a quadratic
+	 * interpolation around its neighbors.
+	 */
+	if (points.size() > 2) {
+		bestPoint = std::clamp(bestPoint, 1ul, points.size() - 2);
+		t = interpolateQuadratic(points[bestPoint - 1],
+					 points[bestPoint],
+					 points[bestPoint + 1]);
+		LOG(Awb, Debug)
+			<< "After quadratic refinement, coarse search has CT "
+			<< t;
+	}
+
+	return t;
+}
+
+void AwbBayes::fineSearch(double &t, double &r, double &b, ipa::Pwl const &prior, const AwbStats &stats) const
+{
+	int spanR = -1;
+	int spanB = -1;
+	double step = t / 10 * kSearchStep * 0.1;
+	int nsteps = 5;
+	ctR_.eval(t, &spanR);
+	ctB_.eval(t, &spanB);
+	double rDiff = ctR_.eval(t + nsteps * step, &spanR) -
+		       ctR_.eval(t - nsteps * step, &spanR);
+	double bDiff = ctB_.eval(t + nsteps * step, &spanB) -
+		       ctB_.eval(t - nsteps * step, &spanB);
+	Pwl::Point transverse({ bDiff, -rDiff });
+	if (transverse.length2() < 1e-6)
+		return;
+	/*
+	 * transverse is a unit vector orthogonal to the b vs. r function
+	 * (pointing outwards with r and b increasing)
+	 */
+	transverse = transverse / transverse.length();
+	double bestLogLikelihood = 0;
+	double bestT = 0;
+	Pwl::Point bestRB;
+	double transverseRange = transverseNeg_ + transversePos_;
+	const int maxNumDeltas = 12;
+
+	/* a transverse step approximately every 0.01 r/b units */
+	int numDeltas = floor(transverseRange * 100 + 0.5) + 1;
+	numDeltas = std::clamp(numDeltas, 3, maxNumDeltas);
+
+	/*
+	 * Step down CT curve. March a bit further if the transverse range is
+	 * large.
+	 */
+	nsteps += numDeltas;
+	for (int i = -nsteps; i <= nsteps; i++) {
+		double tTest = t + i * step;
+		double priorLogLikelihood =
+			prior.eval(prior.domain().clamp(tTest));
+		Pwl::Point rbStart{ { ctR_.eval(tTest, &spanR),
+				      ctB_.eval(tTest, &spanB) } };
+		Pwl::Point samples[maxNumDeltas];
+		int bestPoint = 0;
+
+		/*
+		 * Sample numDeltas points transversely *off* the CT curve
+		 * in the range [-transverseNeg, transversePos].
+		 * The x values of a sample contains the distance and the y
+		 * value contains the corresponding log likelihood.
+		 */
+		double transverseStep = transverseRange / (numDeltas - 1);
+		for (int j = 0; j < numDeltas; j++) {
+			auto &p = samples[j];
+			p.x() = -transverseNeg_ + transverseStep * j;
+			Pwl::Point rbTest = rbStart + transverse * p.x();
+			RGB<double> gains({ 1 / rbTest[0], 1.0, 1 / rbTest[1] });
+			double delta2Sum = stats.computeColourError(gains);
+			p.y() = delta2Sum - priorLogLikelihood;
+
+			if (p.y() < samples[bestPoint].y())
+				bestPoint = j;
+		}
+
+		/*
+		 * We have all samples transversely across the CT curve,
+		 * now let's do a quadratic interpolation for the best result.
+		 */
+		bestPoint = std::clamp(bestPoint, 1, numDeltas - 2);
+		double bestOffset = interpolateQuadratic(samples[bestPoint - 1],
+							 samples[bestPoint],
+							 samples[bestPoint + 1]);
+		Pwl::Point rbTest = rbStart + transverse * bestOffset;
+		RGB<double> gains({ 1 / rbTest[0], 1.0, 1 / rbTest[1] });
+		double delta2Sum = stats.computeColourError(gains);
+		double finalLogLikelihood = delta2Sum - priorLogLikelihood;
+		LOG(Awb, Debug)
+			<< "Fine search t: " << tTest
+			<< " r: " << rbTest[0]
+			<< " b: " << rbTest[1]
+			<< " offset: " << bestOffset
+			<< " likelihood: " << finalLogLikelihood
+			<< (finalLogLikelihood < bestLogLikelihood ? " NEW BEST" : "");
+
+		if (bestT == 0 || finalLogLikelihood < bestLogLikelihood) {
+			bestLogLikelihood = finalLogLikelihood;
+			bestT = tTest;
+			bestRB = rbTest;
+		}
+	}
+
+	t = bestT;
+	r = bestRB[0];
+	b = bestRB[1];
+	LOG(Awb, Debug)
+		<< "Fine search found t " << t << " r " << r << " b " << b;
+}
+
+/**
+ * \brief Find extremum of function
+ * \param[in] a First point
+ * \param[in] b Second point
+ * \param[in] c Third point
+ *
+ * Given 3 points on a curve, find the extremum of the function in that interval
+ * by fitting a quadratic.
+ *
+ * \return The x value of the extremum clamped to the interval [a.x(), c.x()]
+ */
+double AwbBayes::interpolateQuadratic(Pwl::Point const &a, Pwl::Point const &b,
+				      Pwl::Point const &c) const
+{
+	const double eps = 1e-3;
+	Pwl::Point ca = c - a;
+	Pwl::Point ba = b - a;
+	double denominator = 2 * (ba.y() * ca.x() - ca.y() * ba.x());
+	if (std::abs(denominator) > eps) {
+		double numerator = ba.y() * ca.x() * ca.x() - ca.y() * ba.x() * ba.x();
+		double result = numerator / denominator + a.x();
+		return std::max(a.x(), std::min(c.x(), result));
+	}
+	/* has degenerated to straight line segment */
+	return a.y() < c.y() - eps ? a.x() : (c.y() < a.y() - eps ? c.x() : b.x());
+}
+
+} /* namespace ipa */
+
+} /* namespace libcamera */
diff --git a/src/ipa/libipa/awb_bayes.h b/src/ipa/libipa/awb_bayes.h
new file mode 100644
index 000000000000..5284f0ca4cc0
--- /dev/null
+++ b/src/ipa/libipa/awb_bayes.h
@@ -0,0 +1,67 @@
+/* SPDX-License-Identifier: LGPL-2.1-or-later */
+/*
+ * Copyright (C) 2024 Ideas on Board Oy
+ *
+ * Base class for bayes AWB algorithms
+ */
+
+#pragma once
+
+#include <map>
+#include <memory>
+#include <tuple>
+#include <vector>
+
+#include <libcamera/base/utils.h>
+
+#include <libcamera/control_ids.h>
+#include <libcamera/controls.h>
+
+#include "libcamera/internal/yaml_parser.h"
+
+#include "awb.h"
+#include "interpolator.h"
+#include "pwl.h"
+#include "vector.h"
+
+namespace libcamera {
+
+namespace ipa {
+
+class AwbBayes : public AwbAlgorithm
+{
+public:
+	AwbBayes() = default;
+
+	int init(const YamlObject &tuningData) override;
+	AwbResult calculateAwb(const AwbStats &stats, int lux) override;
+	RGB<double> gainsFromColourTemperature(double temperatureK) override;
+	void handleControls(const ControlList &controls) override;
+
+private:
+	int readPriors(const YamlObject &tuningData);
+
+	void fineSearch(double &t, double &r, double &b, ipa::Pwl const &prior,
+			const AwbStats &stats) const;
+	double coarseSearch(const ipa::Pwl &prior, const AwbStats &stats) const;
+	double interpolateQuadratic(ipa::Pwl::Point const &a,
+				    ipa::Pwl::Point const &b,
+				    ipa::Pwl::Point const &c) const;
+
+	Interpolator<Pwl> priors_;
+	Interpolator<Vector<double, 2>> colourGainCurve_;
+
+	ipa::Pwl ctR_;
+	ipa::Pwl ctB_;
+	ipa::Pwl ctRInverse_;
+	ipa::Pwl ctBInverse_;
+
+	double transversePos_;
+	double transverseNeg_;
+
+	ModeConfig *currentMode_ = nullptr;
+};
+
+} /* namespace ipa */
+
+} /* namespace libcamera */
diff --git a/src/ipa/libipa/meson.build b/src/ipa/libipa/meson.build
index c550a6eb45b6..b519c8146d7c 100644
--- a/src/ipa/libipa/meson.build
+++ b/src/ipa/libipa/meson.build
@@ -3,6 +3,7 @@
 libipa_headers = files([
     'agc_mean_luminance.h',
     'algorithm.h',
+    'awb_bayes.h',
     'awb_grey.h',
     'awb.h',
     'camera_sensor_helper.h',
@@ -22,6 +23,7 @@ libipa_headers = files([
 libipa_sources = files([
     'agc_mean_luminance.cpp',
     'algorithm.cpp',
+    'awb_bayes.cpp',
     'awb_grey.cpp',
     'awb.cpp',
     'camera_sensor_helper.cpp',