matrix.hpp

Go to the documentation of this file.

#pragma once
 
// Small matrix operations for spatial audio and graphics.
// Supports 2x2, 3x3, and 4x4 float matrices with common operations.
 
#include <array>
#include <cmath>
#include <cstring>
 
namespace pulp::signal {

template<int N, typename SampleType = float>
struct Matrix {
    std::array<SampleType, N * N> data{};
 
    SampleType& at(int row, int col) { return data[row * N + col]; }
    SampleType at(int row, int col) const { return data[row * N + col]; }
 
    SampleType& operator()(int row, int col) { return at(row, col); }
    SampleType operator()(int row, int col) const { return at(row, col); }

    static Matrix identity() {
        Matrix m{};
        for (int i = 0; i < N; ++i)
            m.at(i, i) = SampleType{1.0f};
        return m;
    }

    static Matrix zero() { return Matrix{}; }

    Matrix operator*(const Matrix& other) const {
        Matrix result{};
        for (int i = 0; i < N; ++i)
            for (int j = 0; j < N; ++j)
                for (int k = 0; k < N; ++k)
                    result.at(i, j) += at(i, k) * other.at(k, j);
        return result;
    }

    Matrix operator*(SampleType scalar) const {
        Matrix result{};
        for (int i = 0; i < N * N; ++i)
            result.data[i] = data[i] * scalar;
        return result;
    }

    Matrix operator+(const Matrix& other) const {
        Matrix result{};
        for (int i = 0; i < N * N; ++i)
            result.data[i] = data[i] + other.data[i];
        return result;
    }

    Matrix transposed() const {
        Matrix result{};
        for (int i = 0; i < N; ++i)
            for (int j = 0; j < N; ++j)
                result.at(j, i) = at(i, j);
        return result;
    }
 
    bool operator==(const Matrix& other) const { return data == other.data; }
};
 
using Matrix2 = Matrix<2>;
using Matrix3 = Matrix<3>;
using Matrix4 = Matrix<4>;
using Matrix2d = Matrix<2, double>;
using Matrix3d = Matrix<3, double>;
using Matrix4d = Matrix<4, double>;

template <typename SampleType>
inline SampleType determinant(const Matrix<2, SampleType>& m) {
    return m(0, 0) * m(1, 1) - m(0, 1) * m(1, 0);
}

template <typename SampleType>
inline SampleType determinant(const Matrix<3, SampleType>& m) {
    return m(0, 0) * (m(1, 1) * m(2, 2) - m(1, 2) * m(2, 1))
         - m(0, 1) * (m(1, 0) * m(2, 2) - m(1, 2) * m(2, 0))
         + m(0, 2) * (m(1, 0) * m(2, 1) - m(1, 1) * m(2, 0));
}

inline Matrix4 translation_matrix(float x, float y, float z) {
    auto m = Matrix4::identity();
    m(0, 3) = x;
    m(1, 3) = y;
    m(2, 3) = z;
    return m;
}

inline Matrix4 scale_matrix(float sx, float sy, float sz) {
    auto m = Matrix4::identity();
    m(0, 0) = sx;
    m(1, 1) = sy;
    m(2, 2) = sz;
    return m;
}

inline Matrix4 rotation_x(float radians) {
    auto m = Matrix4::identity();
    float c = std::cos(radians), s = std::sin(radians);
    m(1, 1) = c;  m(1, 2) = -s;
    m(2, 1) = s;  m(2, 2) = c;
    return m;
}

inline Matrix4 rotation_y(float radians) {
    auto m = Matrix4::identity();
    float c = std::cos(radians), s = std::sin(radians);
    m(0, 0) = c;  m(0, 2) = s;
    m(2, 0) = -s; m(2, 2) = c;
    return m;
}

inline Matrix4 rotation_z(float radians) {
    auto m = Matrix4::identity();
    float c = std::cos(radians), s = std::sin(radians);
    m(0, 0) = c;  m(0, 1) = -s;
    m(1, 0) = s;  m(1, 1) = c;
    return m;
}

struct Vec3 {
    float x = 0, y = 0, z = 0;
};
 
inline Vec3 transform(const Matrix4& m, const Vec3& v) {
    return {
        m(0, 0) * v.x + m(0, 1) * v.y + m(0, 2) * v.z + m(0, 3),
        m(1, 0) * v.x + m(1, 1) * v.y + m(1, 2) * v.z + m(1, 3),
        m(2, 0) * v.x + m(2, 1) * v.y + m(2, 2) * v.z + m(2, 3)
    };
}
 
}  // namespace pulp::signal