Bode.hpp

Go to the documentation of this file.

#pragma once
 
//std
#include <numbers>
#include <cassert>
#include <concepts>
#include <complex>
#include <algorithm>
#include <variant>
 
// eigen
#include <Eigen/Dense>
 
// csvd a csv reader
#include <csvd/csvd.hpp>
 
// controlpp
#include <controlpp/ContinuousTransferFunction.hpp>
#include <controlpp/TimeSeries.hpp>
#include <controlpp/algorithm.hpp>
#include <controlpp/conversion.hpp>
 
namespace controlpp{
 
    template<class T = double>
    class Bode{
        private:
            Eigen::Vector<T, Eigen::Dynamic> omegas_; // frequencies in rad
            Eigen::Vector<std::complex<T>, Eigen::Dynamic> values_; // complex magnitudes
    
        public:
            Bode() = default;
 
            Bode(const Eigen::Vector<T, Eigen::Dynamic>& freqs_rad, const Eigen::Vector<std::complex<T>, Eigen::Dynamic>& values)
                : omegas_(freqs_rad)
                , values_(values)
            {
                // assert same size
                assert(this->omegas_.size() == this->values_.size());
 
                // assert only positive frequencies
                for(size_t i = 0; i < static_cast<size_t>(omegas_.size()); ++i){
                    assert(omegas_(i) >= T(0));
                }
 
                // assert frequencies are ascending
                for(size_t i = 1; i < static_cast<size_t>(omegas_.size()); ++i){
                    assert(omegas_(i) > omegas_(i-1));
                }
            }
 
            Bode(Eigen::Vector<T, Eigen::Dynamic>&& freqs_rad, Eigen::Vector<std::complex<T>, Eigen::Dynamic>&& values)
                : omegas_(std::move(freqs_rad))
                , values_(std::move(values))
            {
                // assert same size
                assert(this->omegas_.size() == this->values_.size());
                
                // assert only positive frequencies
                for(size_t i = 0; i < static_cast<size_t>(omegas_.size()); ++i){
                    assert(omegas_(i) >= T(0));
                }
 
                // assert frequencies are ascending
                for(size_t i = 1; i < static_cast<size_t>(omegas_.size()); ++i){
                    assert(omegas_(i) > omegas_(i-1));
                }
            }
 
            Bode(const Eigen::Vector<T, Eigen::Dynamic>& freqs_rad, Eigen::Vector<std::complex<T>, Eigen::Dynamic>&& values)
                : omegas_(freqs_rad)
                , values_(std::move(values))
            {
                // assert same size
                assert(this->omegas_.size() == this->values_.size());
 
                // assert only positive frequencies
                for(size_t i = 0; i < static_cast<size_t>(omegas_.size()); ++i){
                    assert(omegas_(i) >= T(0));
                }
 
                // assert frequencies are ascending
                for(size_t i = 1; i < static_cast<size_t>(omegas_.size()); ++i){
                    assert(omegas_(i) > omegas_(i-1));
                }
            }
 
            Bode(Eigen::Vector<T, Eigen::Dynamic>&& freqs_rad, const Eigen::Vector<std::complex<T>, Eigen::Dynamic>& values)
                : omegas_(std::move(freqs_rad))
                , values_(values)
            {
                // assert same size
                assert(this->omegas_.size() == this->values_.size());
 
                // assert only positive frequencies
                for(size_t i = 0; i < static_cast<size_t>(omegas_.size()); ++i){
                    assert(omegas_(i) >= T(0));
                }
 
                // assert frequencies are ascending
                for(size_t i = 1; i < static_cast<size_t>(omegas_.size()); ++i){
                    assert(omegas_(i) > omegas_(i-1));
                }
            }
 
            const Eigen::Vector<T, Eigen::Dynamic>& frequencies() const {
                return this->omegas_;
            }
 
            Eigen::Vector<T, Eigen::Dynamic>& frequencies() {
                return this->omegas_;
            }
 
            T& frequency(std::size_t n) {
                return this->omegas_(n);
            }
 
            const T& frequency(std::size_t n) const {
                return this->omegas_(n);
            }
 
            void prewarp_tustin(const T& Ts){
                this->omegas_ = controlpp::prewarp_tustin(this->omegas_, Ts);
            }
 
            void unwarp_tustin(const T& Ts){
                this->omegas_ = controlpp::unwarp_tustin(this->omegas_, Ts);
            }
 
            const Eigen::Vector<std::complex<T>, Eigen::Dynamic>& values() const {
                return this->values_;
            }
 
            Eigen::Vector<std::complex<T>, Eigen::Dynamic>& values() {
                return this->values_;
            }
 
            const std::complex<T>& value(std::size_t n) const {
                return this->values_(n);
            }
 
            std::complex<T>& value(std::size_t n) {
                return this->values_(n);
            }
 
            size_t size() const {return this->omegas_.size();}
 
            bool empty() const {return this->size() == 0;}
 
 
            std::complex<T> value_at(const T& frequency) const {
                std::optional<std::pair<const T*, const T*>> result = find_enclosing(this->omegas_, frequency);
                if(result.has_value()){
                    // extract result values (iterators)
                    const T* low = result.value().first;
                    const T* high = result.value().second;
 
                    // calculate integral ierators
                    const size_t i_low = low - this->omegas_.data();
                    const size_t i_high = high - this->omegas_.data();
 
                    // get frequencies and values
                    const T f_low = *low;
                    const T f_high = *high;
                    const std::complex<T> v_low = this->values_(i_low);
                    const std::complex<T> v_high = this->values_(i_high);
 
                    // calculate percentage for interpolation
                    const T p = (frequency - f_low) / (f_high - f_low);
 
                    // linear magnitude interpolation
                    const T mag_interp = std::abs(v_low) * (static_cast<T>(1) - p) + (p) * std::abs(v_high);
 
                    // linear phase interpolation
                    const T phase_interp = std::arg(v_low) * (static_cast<T>(1) - p) + (p) * std::arg(v_high);
 
                    // calculate resulting complex value
                    const std::complex<T> result = std::polar(mag_interp, phase_interp);
                    return result;
                }else{
                    // assume same value for out of bound values | do not extrapolate
                    if(frequency <= this->omegas_(0)){
                        return this->values_(0);
                    }else{
                        return this->values_(this->values_.size()-1);
                    }
                }
            }
 
            T phase_at(const T& frequency) const {
                const std::complex<T> value = this->value_at(frequency);
                const T phase = std::arg(value);
                return phase;
            }
 
            T phase_deg_at(const T& frequency) const {
                const T phase_rad = this->phase_at(frequency);
                const T phase_deg = phase_rad * static_cast<T>(180) / std::numbers::pi_v<T>;
                return phase_deg;
            }
 
            T magnitude_at(const T& frequency) const {
                const std::complex<T> value = this->value_at(frequency);
                const T mag = std::abs(value);
                return mag;
            }
 
            T magnitude_dB_at(const T& frequency) const {
                const T mag = this->magnitude_at(frequency);
                const T mag_dB = static_cast<T>(20) * std::log10(mag);
                return mag_dB;
            }
    };
 
    template<class T>
    Bode<T> prewarp_tustin(const Bode<T>& bode, const T& Ts){
        Bode<T> result(controlpp::prewarp_tustin(bode.frequencies(), Ts), bode.values());
        return result;
    }
 
    template<class T>
    Bode<T> unwarp_tustin(const Bode<T>& bode, const T& Ts){
        Bode<T> result(controlpp::unwarp_tustin(bode.frequencies(), Ts), bode.values());
        return result;
    }
 
    template<class T>
    const Eigen::Vector<T, Eigen::Dynamic>& frequencies(const Bode<T>& bode) {
        return bode.frequencies();
    }
 
    template<class T>
    Eigen::Vector<T, Eigen::Dynamic> frequencies_hz(const Bode<T>& bode) {
        return bode.frequencies() * std::numbers::inv_pi_v<T> / 2;
    }
 
    template<class T>
    Eigen::Vector<T, Eigen::Dynamic> real(const Bode<T>& bode){
        return bode.values().real();
    }
 
    template<class T>
    Eigen::Vector<T, Eigen::Dynamic> imag(const Bode<T>& bode){
        return bode.values().imag();
    }
 
    template<class T>
    const Eigen::Vector<T, Eigen::Dynamic>& values(const Bode<T>& bode) {
        return bode.values();
    }
 
    template<class T>
    Eigen::Vector<T, Eigen::Dynamic> magnitudes(const Bode<T>& bode) {
        return bode.values().array().abs();
    }
 
    template<class T>
    Eigen::Vector<T, Eigen::Dynamic> magnitudes_dB(const Bode<T>& bode) {
        return static_cast<T>(20) * bode.values().array().abs().log10();
    }
 
    template<class T>
    Eigen::Vector<T, Eigen::Dynamic> phases(const Bode<T>& bode) {
        Eigen::Vector<T, Eigen::Dynamic> result = bode.values().array().arg();
        return unwrap_rad(result);
    }
 
    template<class T>
    Eigen::Vector<T, Eigen::Dynamic> phases_deg(const Bode<T>& bode) {
        Eigen::Vector<T, Eigen::Dynamic> result = bode.values().array().arg() * static_cast<T>(180) / std::numbers::pi_v<T>;
        return unwrap_deg(result);
    }
 
    template<class T>
    TimeSeries<T> impulse(const Bode<T>& bode, const T& time_step, const T& simulation_time){
        assert(bode.empty() == false);
        assert(simulation_time > T(0));
        assert(time_step > T(0));
 
        const T pi = std::numbers::pi_v<T>;
        const std::complex<T> j(0, 1);
 
        const T start_time = 0.0;
        const size_t number_of_samples = static_cast<size_t>(simulation_time / time_step + T(0.5));
        Eigen::Vector<T, Eigen::Dynamic> times = Eigen::Vector<T, Eigen::Dynamic>::LinSpaced(number_of_samples, start_time, simulation_time);
        Eigen::Vector<T, Eigen::Dynamic> values(number_of_samples);
        values.setZero();
 
        const auto f1 = bode.frequencies().head(bode.frequencies().size() - 1);
        const auto f2 = bode.frequencies().tail(bode.frequencies().size() - 1);
 
        Eigen::Vector<T, Eigen::Dynamic> delta_f = f2.array() - f1.array();
        T df_max = delta_f.maxCoeff();
 
        const auto X1 = bode.values().head(bode.values().size() - 1);
        const auto X2 = bode.values().tail(bode.values().size() - 1);
 
        // estimate where to switch from the small t or t=0 solution to the large t solution
        const T t_switch = 0.001 / (2 * pi * df_max);
 
        // complex phase of the current iteration
        Eigen::Vector<std::complex<T>, Eigen::Dynamic> e_j_2_pi_f_ti(bode.frequencies().size());
        e_j_2_pi_f_ti.setOnes();
 
        // complex phase turner (is also the value of the first iteration)
        Eigen::Vector<std::complex<T>, Eigen::Dynamic> e_j_2_pi_f_t1 = (j * T(2) * pi * time_step * bode.frequencies().array()).exp();
 
        const bool includes_dc = bode.frequencies()[0] == 0;
        const std::complex<T> dc_gain = bode.values()[0];
        const T delta_f_dc = bode.frequencies()[0];
 
        // special case t=0
        std::complex<T> v0_one_sided = 0;
        {
            // the one sided result for positive frequencies
            v0_one_sided = ((X1.array() + X2.array()) * delta_f.array()).sum() * T(0.5);
 
            // add artificial sample at 0Hz
            if(includes_dc == false){
                v0_one_sided += ((dc_gain + X1(0)) * delta_f_dc) * T(0.5);
            }
 
            values(0) = std::real(v0_one_sided) * 2; // assume mirrored conjugated negative frequencies
        }
        
 
        // linear approximation (for numerical stability)
        size_t t_itr = 1;
        std::complex<T> LinF = j * (pi/T(3)) * (delta_f.array() * (X1.array() * (T(2) * f1.array() + f2.array()) + X2.array() * (T(2) * f2.array() + f1.array()))).sum();
        
        // add artificial sample at 0Hz
        if(includes_dc == false){
            LinF += j * (pi/T(3)) * (delta_f_dc * (dc_gain * f1(0) + X1(0) * (T(2) * f1(0))));
        }
 
        for(; (t_itr < static_cast<size_t>(number_of_samples)) && (times(t_itr) <= t_switch); ++ t_itr){
            const T t = times(t_itr);
            std::complex<T> v = v0_one_sided + LinF * t;
            values(t_itr) = std::real(v) * 2; // assume mirrored conjugated negative frequencies (imaginary part cancels, real part adds twice)
            e_j_2_pi_f_ti.array() *= e_j_2_pi_f_t1.array();
        }
 
        // exact integration
        for(; t_itr < static_cast<size_t>(number_of_samples); ++ t_itr){
            const T t = times(t_itr);
            const T w = 2 * pi * t;
            
            e_j_2_pi_f_ti.array() *= e_j_2_pi_f_t1.array();
 
            const auto E_f1 = e_j_2_pi_f_ti.head(e_j_2_pi_f_ti.size()-1);
            const auto E_f2 = e_j_2_pi_f_ti.tail(e_j_2_pi_f_ti.size()-1);
            
            const Eigen::Vector<std::complex<T>, Eigen::Dynamic> Ia = E_f1.array() * (T(1) + j * w * delta_f.array()) - E_f2.array();
            const Eigen::Vector<std::complex<T>, Eigen::Dynamic> Ib = E_f2.array() * (T(1) - j * w * delta_f.array()) - E_f1.array();
            const Eigen::Vector<std::complex<T>, Eigen::Dynamic> I = (X1.array() * Ia.array() + X2.array() * Ib.array()) / delta_f.array();
 
            std::complex<T> sum = I.sum();
 
            // add artificial dc sample
            if(includes_dc == false){
                const std::complex<T> Ia_ = (T(1) + j * w * delta_f_dc) - E_f1(0);
                const std::complex<T> Ib_ = E_f1(0) * (T(1) - j * w * delta_f_dc) - T(1);
                const std::complex<T> I = (dc_gain * Ia_ + X1(0) * Ib_) / delta_f_dc;
                sum += I;
            }
 
            const std::complex<T> v = sum / (w * w);
 
            values(t_itr) = std::real(v) * 2; // assume mirrored conjugated negative frequencies
        }
 
        return TimeSeries<T>(std::move(times), std::move(values));
    }
 
    template<class T>
    TimeSeries<T> impulse(const Bode<T>& bode){
        const T max_freq = bode.frequencies()[bode.frequencies().size()-1];
        const T min_freq = bode.frequencies()[0];
 
        const T time_step = static_cast<T>(1) / (T(2 * 4) * max_freq); // oversample 4 times
        const T total_time = static_cast<T>(1) / (T(4) * min_freq); // only use a quarter of the slowest frequency
 
        return impulse(bode, time_step, total_time);
    }
 
    template<class T>
    void integrate(TimeSeries<T>& out, const TimeSeries<T>& in, const T& v0 = T(0)){
        T prev_value = in.values(0);
        T sum = 0;
        out.resize(in.size());
        out.values(0) = v0;
        for(size_t i = 1; i < in.size(); ++i){
            // first order integration
            const T sum_i = (in.values(i) + prev_value) * (in.times(i) - in.times(i-1)) * T(0.5);
            prev_value = in.values(i);
            sum += sum_i;
            out.values(i) = sum;
        }
    }
 
    template<class T>
    TimeSeries<T> integrate(const TimeSeries<T>& in, const T& v0 = T(0)){
        TimeSeries<T> out;
        return integrate(out, in, v0);
    }
 
    template<class T>
    TimeSeries<T> step(const Bode<T>& bode){
        TimeSeries<T> imp = impulse(bode);
        integrate<T>(imp, imp);
        return imp;
    }
 
    template<class T>
    TimeSeries<T> step(const Bode<T>& bode, const T& time_step, const T& simulation_time){
        TimeSeries<T> imp = impulse(bode, time_step, simulation_time);
        integrate<T>(imp, imp);
        return imp;
    }
 
    // operator +
    //-----------------
 
    template<class T>
    Bode<T> operator+ (const Bode<T>& l, const Bode<T>& r){
        assert(l.size() == r.size());
        assert(l.frequencies() == r.frequencies());
        const Eigen::Vector<std::complex<T>, Eigen::Dynamic> result = l.values().array() + r.values().array();
        return Bode<T>(l.frequencies(), result);
    }
 
    template<class T, std::convertible_to<T> T2>
    Bode<T> operator+ (const Bode<T>& l, const T2& r){
        const Eigen::Vector<std::complex<T>, Eigen::Dynamic> result = l.values().array() + static_cast<std::complex<T>>(l)(r);
        return Bode<T>(l.frequencies(), result);
    }
 
    template<class T, std::convertible_to<T> T2>
    Bode<T> operator+ (const T2& l, const Bode<T>& r){
        const Eigen::Vector<std::complex<T>, Eigen::Dynamic> result = static_cast<std::complex<T>>(l) + r.values().array();
        return Bode<T>(r.frequencies(), result);
    }
 
    template<class T>
    Bode<T> operator+ (const Bode<T>& b){
        const Eigen::Vector<std::complex<T>, Eigen::Dynamic> result = +b.values().array();
        return Bode<T>(b.frequencies(), result);
    }
 
    template<class T, int NumOrder, int DenOrder>
    Bode<T> operator+ (const Bode<T>& l, const ContinuousTransferFunction<T, NumOrder, DenOrder>& r){
        const Eigen::Vector<std::complex<T>, Eigen::Dynamic> r_mags = r.eval_frequencies(l.frequencies());
        const Eigen::Vector<std::complex<T>, Eigen::Dynamic> result = l.values().array() + r_mags.array();
        return Bode<T>(l.frequencies(), result);
    }
 
    template<class T, int NumOrder, int DenOrder>
    Bode<T> operator+ (const ContinuousTransferFunction<T, NumOrder, DenOrder>& l, const Bode<T>& r){
        const Eigen::Vector<std::complex<T>, Eigen::Dynamic> l_mags = l.eval_frequencies(r.frequencies());
        const Eigen::Vector<std::complex<T>, Eigen::Dynamic> result = l_mags.array() + l.values().array();
        return Bode<T>(r.frequencies(), result);
    }
 
    // operator -
    //-----------------
 
    template<class T>
    Bode<T> operator- (const Bode<T>& l, const Bode<T>& r){
        assert(l.size() == r.size());
        assert(l.frequencies() == r.frequencies());
        const Eigen::Vector<std::complex<T>, Eigen::Dynamic> result = l.values().array() - r.values().array();
        return Bode<T>(l.frequencies(), result);
    }
 
    template<class T, std::convertible_to<T> T2>
    Bode<T> operator- (const Bode<T>& l, const T2& r){
        const Eigen::Vector<std::complex<T>, Eigen::Dynamic> result = l.values().array() - static_cast<std::complex<T>>(r);
        return Bode<T>(l.frequencies(), result);
    }
 
    template<class T, std::convertible_to<T> T2>
    Bode<T> operator- (const T2& l, const Bode<T>& r){
        const Eigen::Vector<std::complex<T>, Eigen::Dynamic> result = static_cast<std::complex<T>>(l) - r.values().array();
        return Bode<T>(r.frequencies(), result);
    }
 
    template<class T>
    Bode<T> operator- (const Bode<T>& b){
        const Eigen::Vector<std::complex<T>, Eigen::Dynamic> result = -b.values().array();
        return Bode<T>(b.frequencies(), result);
    }
 
    template<class T, int NumOrder, int DenOrder>
    Bode<T> operator- (const Bode<T>& l, const ContinuousTransferFunction<T, NumOrder, DenOrder>& r){
        const Eigen::Vector<std::complex<T>, Eigen::Dynamic> r_mags = r.eval_frequencies(l.frequencies());
        const Eigen::Vector<std::complex<T>, Eigen::Dynamic> result = l.values().array() - r_mags.array();
        return Bode<T>(l.frequencies(), result);
    }
 
    template<class T, int NumOrder, int DenOrder>
    Bode<T> operator- (const ContinuousTransferFunction<T, NumOrder, DenOrder>& l, const Bode<T>& r){
        const Eigen::Vector<std::complex<T>, Eigen::Dynamic> l_mags = l.eval_frequencies(r.frequencies());
        const Eigen::Vector<std::complex<T>, Eigen::Dynamic> result = l_mags.array() - l.values().array();
        return Bode<T>(r.frequencies(), result);
    }
 
    // operator *
    //-----------------
 
    template<class T>
    Bode<T> operator* (const Bode<T>& l, const Bode<T>& r){
        assert(l.size() == r.size());
        assert(l.frequencies() == r.frequencies());
        const Eigen::Vector<std::complex<T>, Eigen::Dynamic> result = l.values().array() * r.values().array();
        return Bode<T>(l.frequencies(), result);
    }
 
    template<class T, std::convertible_to<T> T2>
    Bode<T> operator* (const Bode<T>& l, const T2& r){
        const Eigen::Vector<std::complex<T>, Eigen::Dynamic> result = l.values().array() * static_cast<std::complex<T>>(r);
        return Bode<T>(l.frequencies(), result);
    }
 
    template<class T, std::convertible_to<T> T2>
    Bode<T> operator* (const T2& l, const Bode<T>& r){
        const Eigen::Vector<std::complex<T>, Eigen::Dynamic> result = static_cast<std::complex<T>>(l) * r.values().array();
        return Bode<T>(r.frequencies(), result);
    }
 
    template<class T, int NumOrder, int DenOrder>
    Bode<T> operator* (const Bode<T>& l, const ContinuousTransferFunction<T, NumOrder, DenOrder>& r){
        const Eigen::Vector<std::complex<T>, Eigen::Dynamic> r_mags = r.eval_frequencies(l.frequencies());
        const Eigen::Vector<std::complex<T>, Eigen::Dynamic> result = l.values().array() * r_mags.array();
        return Bode<T>(l.frequencies(), result);
    }
 
    template<class T, int NumOrder, int DenOrder>
    Bode<T> operator* (const ContinuousTransferFunction<T, NumOrder, DenOrder>& l, const Bode<T>& r){
        const Eigen::Vector<std::complex<T>, Eigen::Dynamic> l_mags = l.eval_frequencies(r.frequencies());
        const Eigen::Vector<std::complex<T>, Eigen::Dynamic> result = l_mags.array() * r.values().array();
        return Bode<T>(r.frequencies(), result);
    }
 
    // operator /
    //-----------------
 
    template<class T>
    Bode<T> operator/ (const Bode<T>& l, const Bode<T>& r){
        assert(l.size() == r.size());
        assert(l.frequencies() == r.frequencies());
        const Eigen::Vector<std::complex<T>, Eigen::Dynamic> result = l.values().array() / r.values().array();
        return Bode<T>(l.frequencies(), result);
    }
 
    template<class T, std::convertible_to<T> T2>
    Bode<T> operator/ (const Bode<T>& l, const T2& r){
        const Eigen::Vector<std::complex<T>, Eigen::Dynamic> result = l.values().array() / std::complex<T>(r);
        return Bode<T>(l.frequencies(), result);
    }
 
    template<class T, std::convertible_to<T> T2>
    Bode<T> operator/ (const T2& l, const Bode<T>& r){
        const Eigen::Vector<std::complex<T>, Eigen::Dynamic> result = std::complex<T>(l) / r.values().array();
        return Bode<T>(r.frequencies(), result);
    }
 
    template<class T, int NumOrder, int DenOrder>
    Bode<T> operator/ (const Bode<T>& l, const ContinuousTransferFunction<T, NumOrder, DenOrder>& r){
        const Eigen::Vector<std::complex<T>, Eigen::Dynamic> r_mags = r.eval_frequencies(l.frequencies().array());
        const Eigen::Vector<std::complex<T>, Eigen::Dynamic> result = l.values().array() / r_mags.array();
        return Bode<T>(l.frequencies(), result);
    }
 
    template<class T, int NumOrder, int DenOrder>
    Bode<T> operator/ (const ContinuousTransferFunction<T, NumOrder, DenOrder>& l, const Bode<T>& r){
        const Eigen::Vector<std::complex<T>, Eigen::Dynamic> result = l.eval_frequencies(r.frequencies().array()) / r.values().array();
        return Bode<T>(r.frequencies(), result);
    }
 
    template<class T, int NumOrder, int DenOrder>
    Bode<T> bode(
            const ContinuousTransferFunction<T, NumOrder, DenOrder>& tf, 
            const Eigen::Vector<T, Eigen::Dynamic>& freqs
    ){
        const Eigen::Vector<std::complex<T>, Eigen::Dynamic> complex_magnitudes = tf.eval_frequencies(freqs);
        const Bode<T> result(freqs, complex_magnitudes);
        return result;
    }
 
    template<class T, int NumOrder, int DenOrder>
    Bode<T> bode(
            const ContinuousTransferFunction<T, NumOrder, DenOrder>& tf, 
            Eigen::Vector<T, Eigen::Dynamic>&& freqs
    ){
        const Eigen::Vector<std::complex<T>, Eigen::Dynamic> complex_magnitudes = tf.eval_frequencies(freqs);
        const Bode<T> result(std::move(freqs), complex_magnitudes);
        return result;
    }
    
    template<class T, int NumOrder, int DenOrder>
    Bode<T> bode_hz(
            const ContinuousTransferFunction<T, NumOrder, DenOrder>& tf, 
            const Eigen::Vector<T, Eigen::Dynamic>& freqs_Hz
    ){
        const Eigen::Vector<std::complex<T>, Eigen::Dynamic> complex_magnitudes = tf.eval_frequencies_hz(freqs_Hz);
        Eigen::Vector<T, Eigen::Dynamic> freqs_rad = freqs_Hz.array() * 2 * std::numbers::pi_v<T>;
        const Bode<T> result(std::move(freqs_rad), complex_magnitudes);
        return result;
    }
 
    template<class T, int NumOrder, int DenOrder>
    Bode<T> bode_hz(
            const ContinuousTransferFunction<T, NumOrder, DenOrder>& tf, 
            Eigen::Vector<T, Eigen::Dynamic>&& freqs_Hz
    ){
        const Eigen::Vector<std::complex<T>, Eigen::Dynamic> complex_magnitudes = tf.eval_frequencies_hz(freqs_Hz);
        freqs_Hz *= 2 * std::numbers::pi_v<T>;
        const Bode<T> result(std::move(freqs_Hz), complex_magnitudes);
        return result;
    }
 
    template<class T, int NumOrder, int DenOrder, std::convertible_to<T> T1, std::convertible_to<T> T2>
    Bode<T> bode_hz(
            const ContinuousTransferFunction<T, NumOrder, DenOrder>& tf, 
            const T1& slowest_freq_Hz, 
            const T2& fastest_freq_Hz, 
            const int samples_per_decade=100
    ){
        const T decades = std::log10(static_cast<T>(fastest_freq_Hz)) - std::log10(static_cast<T>(slowest_freq_Hz));
        const T samples = samples_per_decade * decades;
        const Eigen::Vector<T, Eigen::Dynamic> freqs_Hz = Eigen::Vector<T, Eigen::Dynamic>::LinSpaced(samples, std::log(slowest_freq_Hz), std::log(fastest_freq_Hz)).array().exp();
        return bode_hz<T, NumOrder, DenOrder>(tf, freqs_Hz);
    }
 
    template<class T, int NumOrder, int DenOrder, std::convertible_to<T> T1, std::convertible_to<T> T2>
    Bode<T> bode(
            const ContinuousTransferFunction<T, NumOrder, DenOrder>& tf, 
            const T1& slowest_freq_rad, 
            const T2& fastest_freq_rad, 
            const int samples_per_decade=100
    ){
        const T decades = std::log10(static_cast<T>(fastest_freq_rad)) - std::log10(static_cast<T>(slowest_freq_rad));
        const T samples = samples_per_decade * decades;
 
        const Eigen::Vector<T, Eigen::Dynamic> freqs_rad = Eigen::Vector<T, Eigen::Dynamic>::LinSpaced(samples, std::log(slowest_freq_rad), std::log(fastest_freq_rad)).array().exp();
        return bode<T, NumOrder, DenOrder>(tf, freqs_rad);
    }
 
    template<class T, int NumOrder, int DenOrder>
    Bode<T> bode_hz(
            const ContinuousTransferFunction<T, NumOrder, DenOrder>& tf, 
            const int samples_per_decade=100
    ){
        const auto [slowest_freq_rad, fastest_freq_rad] = slowest_fastest_frequencies(tf);
        const T frequency_from_Hz = slowest_freq_rad / (static_cast<T>(10 * 2) * std::numbers::pi_v<T>);
        const T frequency_to_Hz = fastest_freq_rad * static_cast<T>(10 / 2) / std::numbers::pi_v<T>;
        return bode_hz(tf, frequency_from_Hz, frequency_to_Hz, samples_per_decade);
    }
 
    template<class T, int NumOrder, int DenOrder>
    Bode<T> bode(
            const ContinuousTransferFunction<T, NumOrder, DenOrder>& tf, 
            const int samples_per_decade=100
    ){
        const auto [slowest_freq_rad, fastest_freq_rad] = slowest_fastest_frequencies(tf);
        const T frequency_from_rad = slowest_freq_rad / static_cast<T>(10);
        const T frequency_to_rad = fastest_freq_rad * static_cast<T>(10);
        return bode(tf, frequency_from_rad, frequency_to_rad, samples_per_decade);
    }
 
    template<class T>
    void write_csv (std::ostream& stream, const Bode<T>& bode){
        stream << "Frequencies (Hz), Magnitudes (dB), Phases (deg), Real, Imag" << std::endl;
 
        const int n = bode.size();
        for(int i = 0; i < n; ++i){
            const T omega = bode.frequency(i);
            const T f = to_hz(omega);
 
            const std::complex<T> value = bode.value(i);
            const T mag = to_dB(std::abs(value));
            const T phase = to_deg(std::arg(value));
            const T real = std::real(value);
            const T imag = std::imag(value);
 
            stream << f << ", " << mag << ", " << phase << ", " << real << ", " << imag;
            if(i < n-1) stream << "\n";
        }
    }
 
    enum class EBodeCsvReadError{
        CouldNotFindFrequencyVector, 
        CouldNotFindAmplitudeVectors, 
    };
 
    std::ostream& operator<<(std::ostream& stream, EBodeCsvReadError val);
 
    enum class EFrequencyInterpretation{
        AutoHz,     
        AutoRad,    
        ForceHz,    
        ForceRad,   
    };
 
    enum class EMagnitudeInterpretation{
        Auto,    
        ForceAbs,   
        ForceDB,    
    };
 
    enum class EPhaseInterpretation{
        AutoRad,    
        AutoDeg,    
        ForceRad,   
        ForceDeg    
    };
 
    tl::expected<
        Bode<double>, // success
        std::variant<EBodeCsvReadError, csvd::ReadError>> //error
    read_bode_from_csv(
        std::istream& stream, 
        const csvd::Settings& csv_settings = csvd::Settings(),
        EFrequencyInterpretation freq_interp = EFrequencyInterpretation::AutoHz,
        EMagnitudeInterpretation mag_interp = EMagnitudeInterpretation::Auto,
        EPhaseInterpretation phase_interp = EPhaseInterpretation::AutoDeg
    );
 
}