Making your e-drums neighbor-friendly(ish)

Drums are awesome instruments. I mean, you hit stuff with wooden rods and dope sound comes out. It’s the pinnacle of human engineering. 

They are, however, quite noisy, which makes them unsuitable for urban living spaces that cram more people per square-meter than Bangladesh. Luckily, electronic drumkits are a thing, and even more luckily, they don’t sound terrible anymore. Just plug your earphones, and presto, silent druming. Except… If you live in an apartment (as I do), the constant hitting of pads and pedals propagates vibrations through the floor, that will certainly be heard/felt by your donwstairs and nextdoor neighbors. 

The standard solution for this issue seems to be mounting the entire drumkit onto an MDF platform, and placing it on some insulating material: tennis balls are a common alternative, but some people go all in and use sylomer (an industrial-grade dampening polymer) instead. I wanted to save some money and avoid the hassle of getting a slab of MDF into my otherwise not-so-big-apartment (in particular, because the corona-induced lockdown makes it harder to get one cut to size), so I decided to try a 3D-printing route instead.

I own a Roland TD-1DMK kit with a double bass pedal, so I started off by taking all measurements of both the bass and hihat pedals, the diameter of drums feet poles/tubes, and checking up on the nominal sizes for tennis balls. Threw everything into Fusion 360, and voilà. Below is the assembly for the left bass pedal and hihat  – use the “Explode Model” option in the viewer to see the individual parts:

After some ~24 hours of 3D-printing, the end-result turned out nicely:

On the left, the assembly for the bass drum;  center, the hihat – it is connected to the left bass drum, or it’d tip over. On the right, the cups that replaced the drum kit’s original rubber feet. To actually measure the effectivity of the setup, I used a seismometer app on a phone placed ~50cm from the center of the drumkit, directly on the floor. You can see the results below:

As you see, for the tom/snare hits, the attenuation is very effective: you can barely see the vibrations in the scope. The pedals – by far, the biggest culprits – also get a suprising amount of attenuation:

This is the snapshot of both the bass- and hihat-pedal hits. In black, the dampened version, superimposed on the original signal in orange. You’ll notice that the original transients even clipped. From the image, one sees roughly a 1.5:4.5 ratio between the after/before peaks (on Z), which would boil down to a ~70% worst-case improvement. I’ll admit, this is very unscientific, but I didn’t have access to the raw data of the sensors. In any case, I’m actually quite satisfied with the result. Let’s hope the neighbors are too.

’til next time. 

Type-safe scripting with C++ (and other weird explorations)

The why and the what(-the-f***)

Let’s start this short tale with some background. 

For reasons unclear, I’ve started working on my n-th abandonable side-project. Much detail isn’t necessary here: it’s basically a C++ library for performing simple math operations— averaging, sum, standard deviation, autocorrelation and the like. Each operation is implemented as Functor: instantiate it, call it’s operator(), and done! For a concrete example, consider the SumFunctor below (ignore the inheritance for now), which offers overloads to sum two integers or sum (concatenate) two strings

However, there’s a catch: I also want to be able to invoke these functors from some sort of scripting language. The exact language isn’t relevant (though I have my eye on Tcl); suffices to say that it will, as usual, be dynamically typed. The goal here is to suffer be able to play around with these building blocks somewhat faster, potentially being able to integrate them more easily with other languages.  Ultimately, I want to be able to write something like the following: 

Now, here come the juicy questions:

  • since data types in the script will only be known at runtime (i.e., once the script is being executed), how can an hypothetical interpreter going over dispatch the calls to the correct overloads in the C++ functor? And,
  • how can I make a functor’s operator()s invocable in the most type-safe way possible? 

Ideally, I’d like to have something akin to Qt‘s Q_INVOKABLE tag, but I most certainly don’t want to write a Meta-Object Compiler from scratch, and I’d like to avoid devilish macros as much as possible. So, strap in for some weird templated exploration.

Boundary conditions

To better understand where this is headed, let’s define some constraints and solidify some assumptions.

I want an hypothetical interpreter to transparently deal with available functors, regardless of their implementation specifics. Furthermore, I might want to load/unload functions dynamically, e.g. via some sort of import statement. This can be solved quite classically, by having all functors inherit from the same abstract base class (like the AbstractFunctor in the above snippet), which will provide a cohesive minimal API (OOP haters intensify).

How such interpreter accesses the functors isn’t relevant here: they might be stored in a map, vector, list, etc., from which the right one is picked and called. However, there needs to be a cohesive way to pass arguments to the functor, without knowing either their types nor quantities beforehand.  
To approach this, we can take some cues from Python: when writing custom C extensions, data is passed back and forth using PyObjects, which effectively work as containers for runtime data. This is not distant from Qt’s QVariant, — also a generic runtime data container — that appears quite often in the interface between Qt/C++ and QML/JS.
To avoid rolling our own polymorphic container at first, let’s start off by using C++17’s freshly added std::variant instead.

“But”, I hear you say, “std::variant is not runtime polymorphic; all possible types need to be known at compilation time!”. That’s right. std::variant indeed defines a compile-time sum-type. However, it allows for some limited runtime introspection and provides value semantics, all which come in handy. This will temporarily restrict the data types our functors can deal with, but that’s a constraint we shall lift later. For now, let’s define a simple variant_t,

that packages a value coming from the interpreter. Let’s start off by accepting ints, doubles and strings argument types, which shall be bundled in a vector.

Our AbstractFunctor is thus bound to have a method such as void AbstractFunctor::invoke(const std::vector<variant_t>&), through which it can be, well, invoked. Notice that the functors’ return values are being ignored (our invoke method is void, after all). That’s yet another temporary simplification. We’ll get to that. Eventually.

Verifying a variant vector’s veridical validity

At this point, a first step would be asking the question: given a list of arguments (std::vector<variant_t>) and a certain function signature (operator(A, B, C, ...)), how can we verify that the argument types match said signature? 

This turned out to be an interesting point, since it lies somewhere between compile-time constraints and run-time requirements. And, while I could likely concoct some solution exploiting typeid/type_info, I figured it’d make sense to explore a path using templates. That way, I could offload as much type-checking as possible to the compilation phase. 
So, for the SumFunctor::operator(int,int) above, I’d want to write something as matches<int, int>(arglist), which verifies if arglist contains two ints. We’re thus searching for something of the form

where our typename Args... statement makes use of C++11’s template parameter packs. So, strap in for a brief jump into recursive template meta-programming territory!

Tiny trip through templated territory

Template meta-programming is basically functional programming with a cumbersome syntax. In its recursive form, it follows the principles of mathematical induction quite closely: define a base case (a proof for an n = 0) and an inductive step (proof that n = k also holds for n = k + 1), and voilà, the compiler will do the expansions for you. 

For example: in a language better suited for such task, like Haskell, we could compactly define the recursive sum of a list of integers:

Calling something like sum [1, 2, 3] would essentially expand to 1 + sum [2, 3] -> 1 + 2 + sum [3] -> 1 + 2 + 3, unsurprisingly yielding 6. To write something similar in C++, we end up being a bit more verbose:

While the syntax is unfortunately less readable, the base/inductive case pattern is very recognizable. Calling sum(1, 2, 3) spits out 6 as expected. The added constexprs aren’t mandatory, but allow for compile-time computation when possible, which is very neat.
By expanding this C++ example a bit, the sum function can be made to work with basically any combination of types for which operator+ is defined:

Sure, this introduces C++14’s decltype(auto), adds an overkill noexcept for funsies, and throws some std::move and std::forward into the mix, but hey, you get to recursively add apples and oranges together at compile time*. Tasty**!

* Provided operator+(apple, orange) is defined. YMMV.
** Could’ve been even tastier if written with a fold expression, but let’s leave that aside. 

Ok, back to our main point.

Returning to the route of recursive random ramblings

In light of this recently acquired knowledge (and after banging my head against a wall for some time), a fleshed out bool matches<...>(...) function looks like the following:

First off, there is some syntactical sugar to address here: get_if allows me to check at runtime if an underlying std::variant contains a given type. Next, operator sizeof... returns the length of a template parameter pack — we can thus know how many Types were provided as template arguments. Last, but not least, if constexpr is a C++17 godsend, that enables proper conditional branching at compile-time

In the snippet, matches acts as an entry point for our recursive checkArgumentTypeAt, which, as the name indicates, checks if the element at the N-th position of the arglist has the expected type. The base case and inductive step mechanics are the same as the defined in the previous section. Note, however, that different from our recursive sum example, the First and Second types must be explicitly peeled off from the template parameter pack. This is because the signature of checkArgumentTypeAt is the same in both specializations,  and the compiler needs to be provided with enough information for the SFINAE rules to kick in correctly. If we instead had written

since typename... Types can be empty, a call such as checkArgumentTypeAt<1, int>(arglist) would’ve been ambiguous* — as GCC happily confirms:

* Some other type deduction rules are also lurking in the shadows here. For more, here’s a mandatory Scott Meyers talk.

Deterministically dispatching to designated definitions: dope!

Once the matches function verifies the argument list, I want to be able to dispatch it to a corresponding method. At this point, I figured it’d make sense to define the concept of a Dispatcher, i.e., an object responsible for calling a particular operator() overload on our target functor.

The Dispatcher can thus be templated according to the target functor and the signature of the operator() that will be called:

To please the OOP gods (and allow runtime polymorphism), Dispatcher inherits from an AbstractDispatcher with a minimal API. The matches and checkArgumentTypeAt from before are also present, now as const members functions.

Note that Dispatcher has a member f of type Functor. That allows storing the target functor however it is best suited: reference, pointer or RAII member (which provides support for lambdas!). Right below it is an operator(), which will basically invoke the Functor‘s operator() by correctly dereferencing f. In essence, Dispatcher::operator() is just an indirection to leverage C++’s static dispatch rules, and select the correct target Functor::operator() at compile time.

The missing piece of the puzzle is the elusive dispatch method. There are two problems here: 1) how can the N-th element in the arglist be converted to the correct type and, 2) how can Dispatcher::operator() be invoked with the right number of arguments?

Problem 1 is fairly straightforward. If the matches<...>(...) function returned true, it is known that the N-th element in the list has the same type as the N-th type in the Types... template pack. To extract it, one can leverage std::tuple in a very amusing way: 

This TypeAt expression wraps Types... in a  std::tuple, and uses std::tuple_element to fetch the type of the N-th entry — all at compile time. Now, within the Dispatcher, retrieving and converting the N-th argument in arglist is just a matter of writing std::get<TypeAt<N>>( Neat!

A solution for problem 2 builds on top on problem 1, and our complete dispatch method takes form. Unfortunately, it ends up being way less elegant then desired:

Once again, we combine sizeof... and if constexpr to only evaluate the correct conditional branch at compilation time. Unfortunately, the total amount of supported arguments needs to be added by hand, which is rather disappointing. AFAIK, there’s no way to circumvent this (and the fact that even Qt has done something similar further convinces me of that.) If any bright soul sees a better alternative here, send a pigeon

Remaining remarks

Almost. I promise.

When developing a new functor, I’d like to define which methods will be “invocable” via the dispatcher-shenanigans conjured up in the last sections. The aforementioned mythical (but not forgotten) interpreter should see an invoke method that does all the black magic internally. So, let’s circle back to the first code snippet, and define how our AbstractFunctor could look like: 

Neat little detail: different from classes/structs, template parameter pack expansion rules for functions allow packs to appear anywhere in the list (so long as the function signature provides enough cues as to what-goes-where). This can be exploited to write a little makeInvocable helper method that doesn’t require manually specifying the Functor type.
I can then add it to the constructor of our old friend, SumFunctor:

Let’s finally put everything together:

This yields the very unsurprising — but exciting — output:

If you want to play around with this, check/modify/compile the full source-code at Wandbox (or check the GitHub Gist). Have at it!


While this was a fairly long write-up, it ended up being a fairly shallow exploration in both C++ template madness and in the scripting problem outlined in the introduction. There are still quite a few open problems: how to deal with return values in the functors? How to circumvent the limitations of the introduced variant_t? Can we add support to default arguments? The list goes on. 

While a few of those I’ve got worked out already, others are widely unresolved. In any case, I’ll save all these ramblings for a post to come yet in this decade.

Edit: I absolutely cannot NOT mention ChaiScript, which is everything this post wishes to be when it grows up: a fully developed, type-safe scripting language for C++.  Absolutely fantastic code by Jason Turner & Co. If all this nonsense peaked your interest, take a look at their GitHub!

’til next time!

*Alas, additional alliterations aren’t available anymore. 

Bare-metal? Unit-tests? Why not!

Ok, so I’ll try to keep this brief. 

Some context

Unit testing is a good thing, and you should be doing it (don’t worry: that was just me talking to myself, since I miserably fail to follow my own advice). “Well, why?”, you might ask. Simply put, it’s a nice way of ensuring that your software modules/classes/thingies work in the way they should in isolation. Even though coming up with tests takes time, overall this practice tends reduces development complexity and headaches, forces you to write marginally better APIs and serves as a makeshift documentation for your code. 

Stolen from Reddit.

Good. When we’re talking about unit testing for desktop targets, framework options are plentiful. For the C/C++ world, three options pop in mind fairly quickly, but the list is gigantic. Same goes for most of the popular languages that kids use these days. 

However, in the embedded scene, the scenario is different. Computing power tends to be much more limited, and there are often interactions with external hardware that make fully isolated tests hard or unfeasible. If you’re running a *nix box, you can perhaps still profit from the one of the aforementioned frameworks (or try some tailor-made tool), but if you’re dealing with bare-metal targets, you’re pretty much out of luck. On her book, “Making Embedded Systems”, Elicia White skims over some testing practices for embedded systems, but the gist is that there’s no one-size-fits all approach. This thread pretty much echoes this feeling. 

My 2 cents

Getting to the point. As mentioned in another post, I’ve been working with a colleague on a firmware for a VFD. The firmware quickly racked quite a few software components, and keeping track of what worked and what didn’t was getting way to confusing. 

So, I drew some (*ahem*, a lot of) inspiration from Martin Dørum’s unit-test framework Snow – check it out! -, and started working on a bare-bones version of it for the STM32F10x I was using. The result is L_TEST, my attempt at an minimalist unit-test framework for bare-metal targets. The idea behind L_TEST is allowing the user to quickly define test routines and use basic assertion macros to ensure everything is working as it should. Test cases should look straightforward and clear. I really enjoy Dørum’s take on this, so L_TEST employs pretty much the same syntax as Snow does. Take a look at a simple example:

L_TEST uses printf, so you’ll need to implement the _write syscall. The output is color formatted, and should play nice with bash & friends:

Contrary to Dørum’s Snow, L_TEST is unfortunately not header-only, but the implementation is fairly light: 260 bytes for the aforementioned STM32F10x compiled with the -g flag, as shown in the map file below:

The addition of such functions also allowed the macros to be generally simpler, which in turn also reduces code footprint. Still a win, even though you have another source file to add to your list. Speaking of macros, the list still isn’t very plentiful, but here are a few:

  • L_TEST_MODULE(mname, mdesc): Define a test module, with a name and a string description. Test cases should be defined within its body
  • L_TEST_CASE(desc): Define a test case with a description. Test cases can only be defined withen a L_TEST_MODULE body.
  • L_TEST_ASSERT(val): Asserts that the value is non-zero. If the assertion fails, an error message is printed and the test case is aborted.
  • L_TEST_ASSERTEQ_INT(val, ref): Asserts that val == ref, assuming integer values. If the assertion fails, the expected and current values will be printed, and the test case will be aborted.
  • L_TEST_ASSERTEQ_FLT(val, ref): Asserts that the RPD between val and ref is less than 0.1%. If the assertion fails, the expected and current values will be printed, and the test case will be aborted.
  • L_TEST_ASSERTEQ_STR(val, ref): Executes a strcmp between the two supplied strings. If the assertion fails, the expected string and the supplied string are printed.

You can take a look at the full list here.

And thus, 

Unit-testing functionalities like this are just the tip of the iceberg. There’s a lot more that can be done using dummies, mocks, stubs and so on. L_TEST is not, by any stretch of the imagination, an attempt at a full-fledged test framework. However, it did suit my project’s simple needs, so it might suit someone else’s.  You can take a look at L_TEST on GitHub; have your go with it, and let me if you find (or better, when you find) any bugs/issues. Suggestions are welcome! 

’til next time. 


The VFD Series – Part 1: The ups and downs of SPWM

I’ve given up on trying to post periodically here. Let’s see if this brave act of reverse psychology has some effect on my productivity.


Recently, I’ve started working on firmware for a three-phase frequency inverter. While this is absolutely no technological revolution on the grand scheme of things, it’s certainly new ground for me – and deserves some proper notetaking. So, before we dive into anything, let’s do a quick overview. What’s a frequency inverter? Very broadly, a frequency inverter, or Variable Frequency Drive (VFD), is a device that takes a periodic input signal with a certain frequency and generates a periodic output signal with a different, controllable frequency. 

Practically, we’ll be almost always dealing with sinusoidal (or sinusoidal-ish) input and output signals. The reason for that being that mains power is sinusoidal, and most loads we’re interested in (i.e., induction motors) require some form of rotating magnetic field that can be canonically generated by a superposition of sine waves.  Also, on a typical industrial context, mains power is three-phase AC (figure below, left). On this setup, you have three sinusoidal waves, 120º apart. There are a lot of inherent advantages to this arrangement, including a reduction in wire count and gauge, easier hookup to loads and so on. Three-phase AC is a whole world in and of itself, and since I’m not qualified to give you a tour of it, I’ll leave its further comprehension as an exercise to the reader (ElectroBOOM to the rescue).

Left: Three-phase AC power. Right: Simplified three-phase inverter diagram. (Source:

On the above figure, to the right, we see a simplified representation of a three-phase VFD: first, the input three-phase AC signal is rectified into a constant(-ish) DC potential; then, using some form of switching element (e.g., MOSFETs, IGBTs) an output three-phase signal with the desired frequency is generated by combining the output of the U, V and W legs of the circuit. This, in itself, is already the first challenge to be faced: T1, T2, ... T6 are most efficient when operating in the saturation region – i.e., as hard switches, either fully on, or fully off. So, how can we produce a sinusoidal output signal via elements with a binary behavior? 

Mathemagics: the sideways ZOH

First and foremost, by looking at the above diagram, it is very clear that no two switches on the same leg should be simultaneously active at any time – that would configure a short circuit, raise some magic smoke and probably pop a breaker. With that out of the way, we assume that each leg of the VFD can only be in one of two states: on (i.e., upper switch closed, bottom one open), or off (upper switched open, bottom one closed).  This means each phase’s output can be either tied to the DC bus’ voltage, or to the ground/reference potencial. 

So, we are interested in (re)constructing a signal using nothing but switches. Some quick googlin’ points us to a very commonly used technique for signal reconstruction, called Zero-Order Hold (ZOH). This technique allows us to create a continuous-time signal from a discrete-time signal (e.g., a series of numeric values), by holding constant each of its values for an interval T, as shown in the figure below. In a certain sense, this is akin to a Riemann sum, with the area under each sample acting as an approximation to the area under the original signal at that interval. By setting an appropriate T according to the Nyquist-Shannon sampling theorem, the output b_{zoh}(t) signal will contain the harmonic content of the original r(t). Of course, higher harmonics will be present due to the hard transitions between levels, but these should be filtered off. 

In this definition, the signal b_{zoh}(t) has arbitrary, non-discrete values (i.e., the sampled values of r(kT), k \geq 0). However, as mentioned before, the VFD we’re dealing with only allows each phase to be in two discrete states: fully on or fully off. How to circumvent this? While we can’t control the VFD’s amplitude, we can control how long we keep the signal on or off – that is, each pulse can have its width modulated (PWM). Thinking along the lines of the aforementioned Riemann sum, we can try “tipping” each sample of b_{zoh}(t) on its side. So, assume that on an instant t_0, our signal of interest r(t_0) = f_0 (figure below, left). Our reconstructed signal b_{zoh}(t) will hold the value f_0 during the interval T = t_1 - t_0, which results on an area A_0 = (t_1 - t_0)f_0 = Tf_0

Left: b_{zoh}(t) approximates r(t) in the interval [t_0, t_1]. Center: We compute t_s so that A_0 = A_s, assuming f_0 \approx f_sRight: By plotting the relationship between f_s and t_s, \forall r(t), 0 < r(t) < a, we get the sawtooth carrier wave c(t).

Since our VFD signal b_{pwm}(t) can only be 0 or the maximum amplitude a (figure above, center), we wish to find the instant t_s where the area A_s = (t_s - t_0)a approximately matches A_0. This can be written straightforwardly as

    \begin{align*}  A_s &= A_0 \\ a(t_s - t_0) &= f_0(t_1 - t_0) \\ t_s &=\frac{f_0}{a}(t_1 - t_0) + t_0 \end{align*}

By plugging the above expression for t_s into A_s = (t_s - t_0)a, we get A_s = Tf_s. For an adequate interval T, we can safely assume that f_0 \approx f_s, and thus, A_s \approx A_0. Hooray. Now, by also assuming that 0 < r(t) < a, \forall t, we can write the unsurprising relationship between t_s and f_s, in each interval T, as 

    \begin{align*} f_s = a\frac{(t_s - t_0)}{(t_1 - t_0)} && \forall t,  t_0 \leq t < t_1 \end{align*}

This relationship, copied-and-pasted over multiple T intervals, yields the sawtooth-shaped carrier waveform c(t), as drawn in the figure above, right. It now comes as fairly straightforward that b_{pwm}(t) = a when r(t) > c(t) and b_{pwm}(t) = 0 otherwise. With a bit of wit, we can now write b_{pwm}(t) generically as 

(1)   \begin{equation*} b_{pwm}(t) = a (\text{sign} [ r(t) - c(t) ]) \end{equation*}

Neat and compact. This form is also know as natural PMW (or naturally sampled PWM). And, in case you’re wondering: there’s this obscure thing called uniform PWM, but I won’t be touching that. Ever. 

Chop that sinewave… Julienne or Chiffonade?

By plugging the generic sinusoidal wave below

(2)   \begin{equation*} r(t) = R_0 + R_1 cos(2\pi f_1 + \theta_1) \end{equation*}

in our equation 1 above, we get this neat thing called Sinusoidal Pulse Width Modulation (SPWM). Following the intuition developed in the last section, SPWM generates a waveform with the harmonic content of the desired sine wave, as shown in the figure below. I mean, that PWM signal looks like a jumbled mess, but it has the harmonics we’re looking for, believe me. Low-pass-filtering the signal will reveal the modulated sine wave.

The very attentive may have noticed that, in the above picture, the carrier wave c(t) (in green), is not a sawtooth wave as previously defined, but a triangular wave.  In actuality, if we follow the intuition outlined in the previous section, we’ll notice that any triangle-shaped c(t) produces the same A_0 = A_s area equivalence. In practical applications, however, only three different carrier waveforms are used, which yield three basic PWM schemes:

  • Sawtooth: trailing-edge modulation, or left-aligned PWM (figure below, left
  • Triangular: double-edge modulation, or center-aligned PWM  (figure below, center)
  • Inverse sawtooth: leading-edge modulation, or right-aligned PWM (figure below, right)

PWM generation strategies. All waveforms are bipolar and have a 0.2ms period. Left: Left-aligned PWM (sawtooth carrier). Center: Center-aligned PWM (triangular carrier). Right: Right-aligned PWM (inverse sawtooth carrier).

When I faced this carrier-wave-palette for the first time, my initial question was, “ok, cool. So which one do I pick? Is there any difference?”. As it turns out, there’s always some debate over which one to use, but fairly little quantitative approaches to the issue. From an implementation perspective, a triangular carrier wave is a bit more of a hassle. On an MCU, a sawtooth wave can be tipically spawned by a counter that simply counts up and overflows. A triangular wave, on the other hand, requires said timer to go up, then down again. This means that, for a set carrier wave frequency, you have to either provide such counter with a clock signal twice as fast, or sacrifice one bit in the counter’s resolution. But beyond that, is there any modulation strategy that reduces undesired harmonics in the output signal?

This question is relevant for a couple reasons. First, VFD output isn’t usually filtered before it gets to the load – in fact, the load itself acts as a filter. In the case of an induction motor as load, the coils act as RL-filters, smoothing out the input current. Still, lots of undesired harmonic content in the signal might reduce efficiency, produce heat, and cause vibration and audible noise due to magnetostriction

… but if you judge a fish by its ability to modulate a wave …

Cool. So, at this point, the goal is very clear: evaluate the VFD output from different PWM schemes. But evaluate according to what? Let’s pick two neat ones: harmonic content and total harmonic distortion.

Harmonic Content

We are clearly interested in evaluating which harmonics appear on each PWM scheme. My first thought was to look at the C_n coefficients of the compact trigonometric Fourier series, as defined in equation 3 below. 

(3)   \begin{align*}  f(t) = C_0 + \sum_{n=1}^{\infty} C_n \text{cos}(n\omega_0t + \theta_n) \end{align*}

Unfortunately, things are a bit more hairy than that. If we take a look at the PWM equation 1 above, we see the obvious fact the function is periodic in both r(t) and c(t). To apply the Fourier definition above, we’d need to come up with a closed form for a single period of \text{sign} [ r(t) - c(t) ]), which is clearly algebraical masochism (or suicide). In order to analytically analyze the Fourier expansion for such a function, we need to introduce the Double Fourier Series Method: for any function f(x, y), periodic in both x and y, with a period of 2\pi in both axes* we can write:

(4)   \begin{align*} f(x, y) &= C_{00} + \sum_{n=1}^{+\infty}C_{0n} \text{cos}(ny+\theta_{0n})+ \sum_{m=1}^{+\infty}C_{m0}\text{cos}(mx + \theta_{m0}) \notag \\ &+ \sum_{m=1}^{+\infty}\sum_{n=\pm1}^{\pm\infty}C_{mn}\text{cos}(mx + ny + \theta_{mn}) \end{align*}


Well, first off, for all the math inclined folks out there: sorry, but I’m not touching this with a ten foot pole – I’ll be going down the numerical route. However, if you do wish to check out its analytical expansion for various PWM strategies, check this document. Regardless, looking at this expression does provide us with some neat insight about what we should expect to see: the first term on the right-hand side of equation 4 represents a DC component, while the second and third terms represent the harmonics of both y and x, respectively – these are identical to the one-dimensional Fourier Series in equation 3 above. More interesting, however, is the fourth term in that expression. It expresses the sideband frequencies that are generated as a result of the modulation process. We see that n assumes positive and negative integer values, thus yielding upper- and lower-sideband (USB and LSB) spectra around each main harmonic of the carrier frequency. 

*If we want 1 to fit this criterion, we could write it as b_{pwm}(t) = f(x, y) = \text{sign}[r(x) - c(y)], y = 2\pi f_1 + \theta_1, x = 2\pi f_c + \theta_c, where f_c is the carrier’s frequency and f_1 is the modulated sine’s frequency.

Total Harmonic Distortion

While the Fourier series gives us detailed information on the signal’s harmonics, the Total Harmonic Distortion (THD) factor gives us a handy ratio between the harmonics we care about, and the ones we do not. As mentioned, we are interested in producing a pure sine wave, and as such, we care about only a single fundamental harmonic – everything else may be properly labeled as distortion. Our THD can thus be expressed as 

(5)   \begin{equation*} \text{THD}_F = \frac{\sqrt{\sum_{n=2}^\infty v_n^2}}{v_1} \end{equation*}

where v_n is the amplitude of the n-th harmonic and the in “\text{THD}_F” stands for fundamental. Pure sine waves have a \text{THD}_F = 0, square waves have a \text{THD}_F = 48.3\% (percent is a common representation for THD), and higher factors indicate higher distortion – in our case, meaning that less power is going were it should. 


Naturally, we are interested in performing the aforementioned evaluations on the VFD’s output. To that end, our last missing ingredient is a means of properly combining the signals of the U, V and W legs into a single output signal.  As discussed in the introduction, we are interested in three-phase AC inputs and outputs. Such signals can be easily visualized as a phasor projected onto three base vectors, 120º degrees apart** – each projection representing an individual phase. Well, since a picture is worth a thousand words, let the shamelessly stolen GIF below talk for itself:

Vector and time representation of a three-phase system. Amazeballs GIF taken from switchcraft.orgcheck their article out!

In order to properly compute the magnitude of the rotating equivalent vector above (in black), we need to to represent it in an orthonormal basis (as we see above \{U, V, W\} only span \mathcal{R}^2, and hence do not fulfill the criteria for orthonormality). Let’s thus pick the \{\alpha, \beta\} vectors below as our new basis:

A periodic signal v(t) may be represented both in the three-phase \{U, V, W\} base, or in the orthonormal \{\alpha, \beta\} base.

The choice of \{\alpha, \beta\} is arbitrary (as long as they are orthonormal), but done to simplify upcoming calculations (since \alpha represents the real part of the signal, and \beta, its imaginary part).  Now, with a bit of trigonometry, we can represent a periodic signal v(t) shown above in our new base:

(6)   \begin{equation*} \begin{equation*} v(t) = \begin{bmatrix} v_{\alpha} \\ v_{\beta} \end{bmatrix} = %{ \Large \frac{2}{3} } { \small \begin{bmatrix} 1 & \text{cos}(120<span>^{\circ}</span>) & \text{cos}(-120<span>^{\circ}</span>) \\ 0 & \text{sin}(120<span>^{\circ}</span>) & \text{sin}(-120<span>^{\circ}</span>) \end{bmatrix} } = { \Large \frac{2}{3} } { \small \begin{bmatrix} 1 & -1/2 & -1/2 \\ 0 & \sqrt{3}/2 & \sqrt{3}/2 \end{bmatrix} } \begin{bmatrix} v_U \\ v_V \\ v_W \end{bmatrix}  \end{equation*}

This relationship is know as the Clarke transform, and is frequently used in the analysis of three-phase AC circuits. Let’s now assume that the U, V and W phases are each producing individual SPWM signals as per the equation 1, all 120º degrees apart from each other (figure below, right) – notice that the phases’ outputs were normalized to the [0, 1] range. We can now combine them via our definition 6 above, to produce the actual output signal of our VFD (real part \alpha drawn in the figure below, right):

Left: Individual phases of the VFD combined into the output signal via the Clarke Transformation. Each phase’s PWM signal is normalized, so that an off phase is 0, and a on phase is 1. Right: The real-valued part of the VFD’s output signal (\alpha axis).

It is worth noting how each phase is capable of producing only unipolar signals – i.e., signals ranging from 0V to the V_{DC} voltage of the VFD’s DC bus. Their combined output, however, yields true bipolar output. While this might be slightly counter-intuitive at first, imagine all three phases producing a steady PWM signal with a 50% duty cycle. This “balance point” produces zero output (since V_\alpha = 1*0.5 - 1/2*0.5 - 1/2*0.5 = 0, as per equation 6). From this state, by changing the value of one or more phases, we can produce arbitrary output vectors with magnitudes ranging from -2V_{DC}/3 to 2V_{DC}/3. For a bit more discussion on that topic, check this out.

**In case you’re wondering, this 120º offset between the phases is ultimately related to the physical placement of the stator windings inside three-phase motors and generators.

Last, but not least

We are finally ready to answer the question we posed several paragraphs ago: which PWM strategy is best? Left-aligned, right-aligned or center-aligned? 

First off, as we’ve discussed before, our evaluation tools only care about the frequency spectra of our signals. So, taking into account the time reversibility of the Fourier series, and noting that left- and right-aligned SPWM waveforms are mirror images of each other, we immediately know that they’re equivalent (for our intents and purposes). So, we’ll only compare trailing-edge and double-edge modulations. 

To generate the SPWM signals, I’ve implemented the SPWM definion (i.e., combining equations 1 and 2) in Matlab (Code!). Computing the Fourier single-sided amplitude spectrum (Code!), as well as computing \text{THD}_F (Code!) was also done in Matlab. [Any feedback on the correctness of these code snippets would be greatly appreciated] The generated plots are shown below. U, V, and W phases are modulating a 50Hz sinewave with an amplitude modulation index of 0.8, on a 5kHz carrier (identical as shown in the figure above, left), so the expected output amplitude waveform is 0.8*0.5 = 0.4:

Top and Bottom show Trailing- and Double-Edge Modulation schemes, respectively. Left: Time plots of the SPWM signals. Output is a 50Hz sine wave on a 5kHz carrier. Right: Single-sided power spectra for each signal, with added \text{THD}_F factor.

And, voilà. We immediately see that in both modulation schemes, the desired fundamental is there, almost perfectly in the desired amplitude (0.39 \approx 0.4) – yey! all that SPWM hassle does work after all. Now, interestingly, we have a rather curious result with the FFT spectra and the THD factors. At first, we see that Trailing-Edge modulation has a somewhat smaller distortion factor (70.13%), but its spectrum seems arguably more messy. Moreover, we can see that Double-Edge modulation seems to have much less harmonic content (in fact, it has exactly half of the sideband harmonics, as verifiable if we expand equation 4 analytically for each scenario, as seen here). On top of that, the fundamental switching harmonics (around 5kHz) have smaller amplitudes in the Double-Edge scheme. So, what gives?

It seems that the higher-order harmonics of the Double-Edge modulation tend to weight-in more heavily in the quadratic sum of the THD factor, yielding a higher overall distortion (88.94%). In practice, however, the RL-filter-characteristic of VFD loads will have a cutoff frequency around the hundreds of Hz, so realistically, harmonics above the switching frequency will have almost no effect***. So, we can confidently argue that, in practical applications, Double-Edge modulation – a.k.a. center-aligned PWM – does produce less harmonics, which seems to echo the faint opinions on the topic that float around the interwebs.

Now, of course: as the image above shows, the difference isn’t all that extreme, and as we’ve discussed above, there’s a bit more implementation effort associated with center-aligned PWM. So, once again – and slightly disappointingly – YMMV.

***Very broad and oversimplified generalization. Don’t sue me if you fry your setup testing something I’ve said. 

Disclaimer (& closing thoughts)

Well, let’s just make it very clear: this whole thing is a somewhat brief write-up of my latest incursions in what’s an unknown territory for me. I’ve been figuring stuff out on the fly, so if you spot anything wrong, please, let me know

Edit: In a recent conversation, a friend of mine and a literal master of all-things-electric, Julio, added some very relevant information to this mix. He confirmed that, in practical applications, the choice of carrier waveform is not of great impact. However, when implementing a VFD (e.g. on a MCU) with any kind feedback control loop, the peaks and valleys in the triangular carrier of the center-aligned PWM can be used to synchronize ADC samplings of the generated waveform (figure below). This ensures that the sampling doesn’t happen during switching, reducing measurement noise (and, in case of current measurements, providing that pulse’s average value). This article goes into more depth into that, and it’s worth taking a look. He also mentioned that sawtooth carriers (in trailing- and leading-edge modulations) essentially synchronize the switching in all phases. This increases output noise due do parasitics in the circuit/load (stuff that we did not capture in this write-up), and can be a real issue in high-power applications. Thanks for the insight, Julio! 

Using peaks and valleys of the triangular carrier to synchronize the ADC sampling of the PWM signal, as to avoid switching noise (since measurement happens exactly in the “middle” of an On or Off state). Image source:

’til next time. 

AS5030 magnetic encoder: capturing a PWM signal with an ATSAMD21

It seems that I can’t avoid periodically ostracizing this page. Welp, let’s try to make it up.

Context: the AS5030 magnetic encoder IC

In a project I’m currently working on (more about it in later posts, perhaps), I needed a halfway decent way of measuring the angular displacement of a small, manually-turnt wheel. I had originally mounted a small mechanical encoder (Bourns PEC11R) into my solution, which yielded about 96 pulses per-revolution (PPR). However, on a long run, contacting encoders are not the best pick for these applications due to mechanical wear.

Searching for a better solution I came across ams’ portfolio of magnetic encoders, which features interesting solutions for contactless position measurement. Amongst them, the AS5030 (datasheet) was the one that best met my requirements, providing an improved 256 PPR. The gist of it is simple: mount the IC, power it up and spin a magnet on top of it. The IC will generate a PWM signal with a pulse width proportional to the magnet’s angle in relation to the chip. Spoiler alert: not any magnet will work! You’ll have to get diametric magnets, such as this little guy here.

AS5030 overview: the hall frontend of the chip measures the orientation of the field lines of a diametric magnet and produces both PWM and analog signals accordingly. A SPI interface is also available.

As you may see in the diagram above, the AS5030 also provides a serial interface. The thing, however, is a weird half-duplex SPI that operates on 21-bit long frames, making it odd to use. On top of that, the AS5030 is strictly 5V compatible, which means I would have to level-shift all the signals going to my 3v3 microcontroller. That would just add unnecessary parts and complexity, so I ditched the serial interface and went with PWM (and a resistive divider to shift the voltage for the uC).

The PWM signal frequency is rated at 1.72kHz, but may vary slightly with temperature. The duty cycle encodes the position, going from a spec’d minimum of 2.26us to a spec’d maximum of 578.56us, like shown in the picture below:

AS5030 PWM signal specifictaion.

Getting the angular position accurately can be done via the ratio of the duty cycle (t_{\text{on}}) and the PWM’s period (t_{\text{on}} + t_{\text{off}}), as per the equation supplied in the datasheet:

\text{angle}[{^\circ}] = \frac{360}{256}\large[\large(257\frac{t_{on}}{t_{\text{on}}+t_{\text{off}}}\large)-1\large]

 ATSAMD21: peripheral bureaucracy

Measuring the pulse-width of a PWM signal is a textbook example of input capture. Capture operations allow you to record a signal edge together with a timestamp directly via hardware, without having to employ any CPU cycles. The vast majority of modern microcontrollers support this feature, and the ATSAMD21E18A employed in this project is no exception.

As we see in the ATSAMD21 datasheet (section 31.2), the TCC peripheral not only supports this sort of operation, but also has a dedicated mode for pulse-width capture, where “period T will be captured into CC0 and the pulse-width t_p into CC1″. Unfortunately, there’s a catch. Usually, capture operations occur entirely within a timer peripheral, which detects the signal’s edge on a dedicated pin an stores the timestamp based on its internal counter. However, to increase flexibility, capture operations in the ATSAMD21 use the External Interrupt (EIC) peripheral to generate an event in the internal Event System (EVSYS), which gets relayed via the internal event channels to the Timer/Counter (TCC) peripheral, in turn triggering the capture. Wait, what? Let me sketch that out:

Signal flow on ATSAMD21 capture operation: the PWM signal is routed from the GPIO pin to the EIC peripheral, generating an EVENT (IRQ-esque feature) on high logic. The event edges get relayed via an EVSYS channel, where the EIC acts as a generator, and the TCC timer peripheral is the event user.

This certainly allows for a lot more flexibility, not tying the capture operation to a particular pin, since all pins have EIC Interrupt/Event capability. However, configuring it is quite a mouthful (and quite poorly documented, BTW).

Talk is chip, show me the encode

I suck at puns. Ok, so let’s get a bloatw…I mean, ASF-free setup going. In my particular example, I’ll be using pin PA11, tied to EIC channel 4, but again, any pin can be used. Also, I’ll be using EVSYS’s channel 0, but you can use any of the available channels for any event. This will be done in four steps:

We’ll first configure the timer. It’ll run off the main GCLK (in this case, at 48MHz), counting up from 0 until 0xFFFFFFFF (then wrapping around):

Then, let’s setup the EIC. Sense is set to HIGH, and a tiny helper function configures the channel. Notice how the EVCTRL bit is set, which generates the EIC event:

Let’s then configure the EVSYS: EIC acts as a event generator on channel 0, and TCC1 is the event’s user (akin to a listener):

Last but not least,  let’s configure the pin. The PINMUX_PA11A_EIC_EXTINT11 value should be defined in the samd21.h header:

That’s it! Your device is now running PWM pulse-width capture on pin PA11, with no CPU cycles being used for it at all.

Now what?

Now that everything’s configured, … Well, make sure the PWM signal is getting to the designated pin. The signal’s pulse-width and period are now constantly being fetched and saved into the CC0 and CC1 registers, respectively. To wrap it up with the AS5030, we can now compute the measured angle with the following function:

This function returns angles x10, (i.e. 15.4° = 154), so that no float support needs to be added (which can eat a lot of Flash on smaller devices). It’s also enough to properly deal with the AS5030’s 1.4° resolution.

Best part? It actually works:

Rotating a magnet in front of the AS5030. Very professional test setup.

You can optionaly use the TCC_INTFLAG_MC1 interruption if you don’t want to poll the registers for changes. Yeah.

’til next time.

Getting OpenProject up and running

Ok, so I was looking for a self-hosted project management tool. Something that would fill the gap that the late (and discontinued?) DotProject left in my heart. After lots of Googlin’, I came across OpenProject. Though they offer hosting plans, you can host the tool for free on your own server. Installation looked like it had a bit too many steps for my own taste, but I went with it. Needless to say, if I’m writing this note-to-self, sh*t went south. Considering my overall stupidity regarding all-things-web, I found the installation process fairly convoluted (and apparently I’m not alone). The steps I took to get it working are sumarized below (which were compiled from quite a few different sources).

Rolling up the sleeves

A little context here – I’m installing OpenProject version 6 (latest stable release at the time of writing) on a VM running Ubuntu 14.04. Unfortunately, for other distros, YMMV. For many distros (Ubuntu included), there’s the possibility of going through with the packaged installation (e.g. via apt-get), but I ended up following the manual process. I’ll be referring to these instructions during this post. If the link breaks or changes in the future, get the PDF version (but these instructions will probably be worthless anyway):

Let’s start with the Prepare Your Environment part. I’ve totally skipped the groupadd stuff, and installed everything under my regular user – do as you please. I installed all essentials and Memcached  as listed. I already had a working install of MySQL, so I skiped the Installation of MySQL part, including the database creation (which will happen automatically later on). I then followed the Installation of Ruby (via rbenv)…

…and Installation of Node (via nodenv, but using the 6.0.0 LTS version) as listed on the manual :

Adding some swap and getting to it

I had trouble getting the rbenv install to complete. It would get killed halfway through, in which seemed an insuficcient memory issue. Anticipating a future instruction from the manual, I’ve decided to try adding some swap space to my VM. As usual, DigitalOcean has a great documentation page on how to do this. Spoilers:

Ok, now fetch the relevant version of OpenProject – I’m going with the regular openproject (not the Comunity Edition), version 6, and install it:

Now run the Configure OpenProject and Finish the Installation of OpenProject as described in the instructions. After getting them done, you’ll have an openproject/public folder that looks like this:

The lack of an index file looked super weird to me, but it’s supposed to be that way. Don’t panic.

Getting it to play nice with Apache

I started having some trouble in the Serve OpenProject with Apache section of the instructions. The compilation and installation of the passenger module worked nicely, and I could a2enmod passenger with no trouble. However, the supplied VirtualHost configuration file didn’t cut it for me.

I wanted to have OpenProject accessible from, and since my VM already serves my WordPress and Pydio, I couldn’t just move my DocumentRoot around. So I edited my  configuration in /etc/apache2/sites-enabled/ as follows:

I kept my DocumentRoot where it was, and served OpenProject through an Alias. Setting the PassengerResolveSymliksInDocumentRoot is necessary here since Passenger won’t solve symlinks automatically for versions above 2.2.0 (which happens to be the case). Also, we have to point the PassengerAppRoot to where the app’s is stored — in this case, the root of OpenProject’s cloned git repo.

Also, I’ve added rails_relative_url_root = "/op" to the config/configuration.yml created during the Configure OpenProject step, to match the folder alias I’ve created.

Reload Apache with a sudo service apache2 restart and give it a go. Works on my machine.

’til next time.

openVPS: Poor man’s motion capture

The Why and The What

Flying tiny drones indoors is cool, no questions asked. And stuff gets all the more interesting when you can accurately control the drone’s position in space — enabling all sort of crazy manouvers. However, using regular GPS for such applications is often a no-go: antennas may be too heavy for the drone, accuracy is above the often-desired cm range and getting GPS to work reliably indoors often requires breaking some laws of physics. For that reason, labs and research teams often turn to commercial marker-based motion capture systems. You place fancy cameras around a desired scene, fix some reflective makers on the objects you want to measure/track and some magic software gives you sub-centimeter (often sub-millimeter!) tracking accuracy. Yey. Too bad such systems are expensive. Like, way too expensive. Like, “I’d have to sell my organs on the black market for that”-expensive.

I thought: “surely someone wrote a piece of FOSS that allows me to grab some 3 or 4 webcams I have lying around and setup a poor-man’s version of such a system”. To my amazement, I was wrong (I mean, Google didn’t turn up anything, but I’d still love to run into something).  Considering that I still had to come up with some image-processing-related activity for my Masters degree, this seemed like a nice fit. So, let’s get on with it. [Disclaimer: heavy usage of OpenCV ahead.]

First things first: what are we looking for? Well, the sketch bellow sums it up: we have N cameras (C_0, C_1, … C_{n-1}), all looking at some region of interest from different angles. When an object of interest \{O\}  (with some red markers) shows up in the field of view of 2 \leq k \leq N cameras, we can intersect the rays emanating from each of them to discover the object’s position and orientation in space.

Sketch of system’s layout.

I’m a simple camera – I see the object, I project the object

Ok, but how? Well, let’s understand how a camera works. The image below (which I shamelessly stole) represents a pinhole-camera, a simplified mathematical camera model:

Pinhole camera model.

Any point P = (X,Y,Z) in the scene can be represented as a pixel \alpha = (u,v) via the transformation in equation (1) below. [R|t] are the rotation and translation between the camera’s frame of reference and the scene’s.

(1)   \begin{equation*} s \begin{bmatrix} u \\ v \\ 1 \end{bmatrix} = \begin{bmatrix} f_x & 0 & c_x \\ 0 & f_y & c_y \\ 0 & 0 & 1 \end{bmatrix} \bigg( R \begin{bmatrix} X \\ Y \\ Z \end{bmatrix} + t \bigg)\end{equation*}

The 3×3 matrix on the right-hand side of the equation is called the intrisic parameters matrix, and includes information about the camera’s focal distance and centering on both the and axes. As I’ve mentioned, the pinhole model is a (over)simplified one, which doesn’t accout for distortions introduced by the lens(es). Such distortions can be polynomially modelled, and the coefficients of such polynomial are then addressed as the distortion coefficients. No need to get into detail here, but both intrisic parameters and distortion coefficients can be obtained through OpenCV’s calibration routines. We’ll now assume all cameras are properly calibrated and lens distortion have been compensated, making the pinhole model valid.

Of blobs and vectors

Ok, so first thing: how to find the markers — reflective or not? Well, there’s lots of room to explore and digress here, but I’m lazy and just went with OpenCV’s built-in SimpleBlobDetector (great tutorial here!). Setup and usage are straightforward and results are nice:

Detecting markers (a.k.a. blobs).

Ok, so the detector spits out center coordinates (in pixel space, along the \{u,v\} axes) of each blob’s center point. If we want to “intersect” each view of each marker as previously sketched, then we need to compute unit vectors that emanate from each camera’s optic center \mathcal{F}_c through each blob’s center. Lets assume a point p = (x,y,z) – e.g., the marker – already in camera’s coordinates. Since the projective relationship in equation (1) loses depht information, we’ll rewrite the point in its homogenous form p_h = (x/z,  y/z, 1). We can then compute our unit vector \hat{p} of interest from the blob’s center pixel \alpha_h = (u,v,1) via,

p_h = F^{-1} \alpha_h

\hat{p} = \frac{p_h}{|p_h|}

is our aforementioned intrisic parameters matrix. Implementing this is as straightforward as it looks:

Where art thou, corresponding vector?

Following our system sketch up there, the next logical step is to intersect all these unit vectors we computed in each camera, for each marker. But the question that arises is: how do we know which vectors to match/intersect? — i.e., which vectors are pointing to the same marker on each camera? One convenient way of tackling this is exploring epipolar geometry (thanks for the inspiration, Chung J. et al). Simply put, knowing the rotation and translation between two cameras – e.g., and O’ (see image below) – allows you to map a point X as seen by O to an epiline l’ in O’. This drastically reduces the search space for the corresponding X’ — that is, the corresponding view of the same object.

There’s no need to delve into epipolar geometry — specially since I didn’t learned it in detail. OpenCV supplies satisfactory documentation and great features for this. First, I’ve used OpenCV’s stereoCalibrate within a routine that continuously checks which pairs of cameras see the calibration pattern at any given time, records that information, then later computes the fundamental and essential matrices. These which encode the [R|t] relationships between cameras, as shown in the figure below. Following that, computeCorrespondEpilines get’s you parameters a, b and c for the epiline’s equation au + bv + c = 0. So, if the n-th camera sees marker j, and a camera n+k, k \geq 1 sees a set of markers, the matching view of is m, as given by

\operatorname*{argmin}\limits_{m \in M} au_m+bv_m+c

Encoding this relations yields the results shown below (step-by-step usage of the aforementioned OpenCV functions can be found in the provided stereo_calib.cpp sample – worth checking out).

Epipolar relations between cameras being used to match points between two scenes.

Transformations between cameras

So far we have our unit vectors represented on each camera’s coordinate frame, and we are able to match the vectors “pointing” to the same marker/object on the scene. Great! Now, as a last step before computing their intersection, we need to express them all in a single coordinate frame. For simplicity, the first camera (n = 0) was chosen as the main reference for the system. So, in order to express all unit vectors in that coordinate frame, we need to know [R_n^0|t_n^0],\text{ }\forall n. In theory, we could concatenate all the successive transformations obtained with the stereo calibration procedure, but I’ve foud that it propagates lots of measurement errors and uncertainties. For now, since I’m using a small rig with only 4 cameras, the solution is to simultaneously show all cameras a known reference (namely, an AR marker), as show below:

Computing relationships between cameras via simultaneous imaging of an AR marker.

With each camera seeing the marker at the same time, we can easily compute the relationship between the 0-th camera and the n-th camera by just composing two transformations, as illustrated below:

Transformations between the AR marker, the 0-th camera and the n-th camera.

We can then easily obtain the transformations we’re interested in:

t_n^0 = -R_n^0 t_n^A + t_0^A   R_n^0 = R_0^A(R_n^A)^T

Getting to the point

Now we (finally!) have everything we need to intersect our vectors and find the position of the point(s) in the scene. I greatly recommend you take a look at this post in the Math StackExchange, since I’m piggybacking on Alberto’s detailed answer.

So, let’s assume that k \in [2, N]  cameras see a marker – thus we have unit vectors to deal with. Let’s also assume from now on, that all points and vectors are being represented in the 0-th coordinate frame. Essentially, by “intersecting” the vectors, we are trying to find a point \tilde{P} that is as close as possible to all the lines, minimizing the quadratic error relation given by

(2)   \begin{equation*}  E^2({\tilde{P}}) = \displaystyle \sum_k \bigg( ||(\tilde{P}-\mathcal{F}_{ck})||^2 - \big[ (\tilde{P} - \mathcal{F}_{ck}) \cdot \hat{p_k} \big]^2 \bigg) \end{equation*}

where \mathcal{F}_{ck} is the focal point of the k-th camera and \hat{p}_k is the unit vector emanating from the k-th camera towards our marker of interest. This, ladies and gentlemen, is also probably the most convoluted way of writing the Pythagorean Theorem. To minimize this relationship, we look for its inflection point

\frac{\delta E^2(\tilde{P})}{\delta \tilde{P}} = 0

Taking the derivative of (2) and manipulating it a bit, we get

\displaystyle \bigg[ \sum_k \big[ \hat{p_k} \cdot \hat{p_k}^T - I \big] \bigg] \cdot \tilde{P} = \displaystyle \sum_k \big[ \hat{p_k} \cdot \hat{p_k} \big] \cdot \mathcal{F}_{ck}

This can be represented as the matrix system S \tilde{P} = C. Thus our intersection point is given by \tilde{P} = S^{+} C, where S^{+} = S^T(SS^T)^{-1} is the pseudoinverse of S.

But will it blend?

How disappointing would it be, not to have a single video after all this boring maths? There you go.

If you like looking at badly-written, hastily-put-together code, there is a GitHub page for this mess. I plan on making this a nice, GUI-based thing, but who knows when. For more info, there’s also a report in portuguese, here.

’til next time.

libFilter add-ons

Following my last post on my minimalist filter library, I just got off my butt to add some high-pass filtering capabilities too. That’s specially useful when you’re trying to remove trends from datasets – this happens a lot for instance in biomedical applications, e.g. ECG, where some breathing artifacts come up as low frequency trends on your signal.

A simple example of such trend removal is included. As shown in the snippet below, taken from the example, after instancing the Filter, extracting the trend and DC components is very straigthforward:

The result looks like the following (blue senoidal wave on top is signal + trend + DC, red senoid is signal + DC, orange senoid is the filtered signal):

The filters implemented go only up to order 2 (I was lazy, sorry), and are also based on normalized Butterworth polynomials, discretized via Bilinear Transforms. Frequency responses for two hypotetical filters at 1kHz sampling frequency are shown below:

All filters are implemented in closed forms, so they can be reparameterized on the fly with minimum computational effort. Enjoy (and report any bugs). Major cleanup on the code is foreseen (and a notch-filtering example too).

’til next time.

Minimalist low-pass filter library

So, the other day I needed to compute some low-pass filters on the fly on an Arduino Mega. By “on the fly” I mean that the filters’ parameters would eventually be recomputed mid-operation, so setting some equation with static const parameters would not cut it. Using a static filter, is, however, the most common application scenario. If that’s your case – and I insist – tweak the sh*t out of your filter with a decent tool then just add it to your code. If, on the other hand, you need to update your filter dynamically (or if you’re just plain lazy to compute poles and zeros), then this is for you.

I ended up putting together a minimalist library for that, libFilter (github!). Works with the Arduino IDE and as a regular C++ library (just uncommenting one #define does the trick).

Using two filter instances on a signal from a load cell getting punched.

Using two filter instances on a signal from a load cell getting punched.

For now it only implements low-pass filters based on  normalized butterworth polynomials, but who knows what necessity might add to it next. Let’s take a look at some ADC filtering:

’til next time!

P.S.: I just can’t avoid letting this page eventually fall into temporary ostracism. Geez.

Finding my way with Cura 10.06

I’ve been using Cura as my go-to 3D-printer slicer for quite some time now. Compared to Slic3r, it’s faster and produces more optimized G-Code (using the same settings in both slicers, Cura’s prints were faster for me- but as always, YMMV). However, Cura provides less tweakable options than Slic3r, so it takes some getting used to.

So, I was using Cura 10.04 under OSX (10.10 and 10.11), and all was nice and good. Suddenly, I wild bug appeared. As shown in the picture below, Cura uses some text overlays on the 3D-viewer to display information: print duration, filament usage, estimated cost, type of view, etc. On my machine, all text overlays were suddenly gone! Rebooting, reinstalling, re*ing, etc, didn’t seemed to solve the issue. And, needless to say, using a slicer without the aforementioned infos is quite frustrating.


Left: Fully functional Cura, with text overlays in 3D-view (image courtesy Tom’s Guides). Right: My wild bug in action – bye, text overlays.

Furthermore, this bug seemed to afflict only a handful of unlucky bastards: Google pointed me to one unanswered forum entry, and to one GitHub issue ticket. The latter got a response from one of Cura’s developers, which stated that they’re “not working on a fix, as we will be repacing the whole UI in a few months, which no longer uses this old font rendering code”.

All right then. The dev’s answer was posted mid-April 2015, so I figured a new version should already exist. Indeed, in the list of beta releases, Cura 10.06 was available for download. I grabbed it and got going.

The UI got heavily revamped, and lots of options were made available (bringing a Slic3r-y level of tweakability to Cura). A quick tour in Cura’s Preferences > Settings allows you to select which options you want to edit in the main window. No need for much detail here, as the help pop-ups on each option make the whole experience nicely intuitive.


Cura 10.06 revamped UI.

Unfortunately, as new options were added, others were lost. When first configuring Cura, you need to add a printer. In previous versions, you could choose from a list of models or just go with a “generic” machine – to which you’d add dimensions and other relevant informations. This option is no longer readily available – you have to pick a model from the list! Having a DIY printer that looks like nothing on the list, I was quite frustrated. For test purposes, I picked the Prusa-i3 preset (though my printer uses a bed roughly four times as big).

On top of that – and that was the killer for me – there is no option to edit the Start/End G-Code procedures. While this sounds trivial/useless for most off-the-shelf printers, it won’t (probably) work with a custom machine (e.g., like mine). For instance, due to the form and size of my printer’s bed, I can’t go homing XYZ like a Prusa-i3 would do (and as the preset G-Code insists on doing).

After shelving Cura 10.06 for some days, I stumbled on this page on the RepRap Wiki. It shed some precious light on the workarounds required to solve these (and other) problems. As it turns out, much of Cura’s printer configurations are set in JSON files – a file per printer model, stored in Cura’s installation folder, under resources > machines. To add my custom printer, I’ve copied the prusa_i3.json preset to a new custom_3dcnc.json in the same folder, and went on editing it. The JSON entries are pretty self-explanatory:

Changing the id and name fields makes your custom printer appear on Cura’s list. Also note: AFAIK, you need a STL model for your printer’s bed (commenting out the field platforms won’t work). By default, they’re stored under resources > meshes. I could reuse one, but I’ve simply exported my printbed’s STL from my CAD program.

Configured Cura (with my printer's crazy-ass printbed on display).

Configured Cura (with my printer’s printbed being displayed).

Last, but not least, machine_start_gcode and machine_end_gcode are what I was looking for. Add the G-Code inside a single pair of quotes, with commands spaced with a newline character (\n) and you should be golden. Save the file, reload Cura, and you’re good.

’til next time.