Orders of complimentary butterworth filters for IMU I'm currently trying to design a complimentary filter for deployment on a drone, as part of an assessed semester project. Even though processing power is not an issue (we have GHz, lots of pipelining and all the good stuff), I am still electing to use an IIR filter in the final implementation, as this is easier to represent mathematically in my opinion.

As I have read Here in this answer, certain orders of filters do not produce monotonic unity gain in output. As such, at this stage I am mostly interested in proving, in a round-about way, which orders of filters produce monotonic unity gain. Phase difference is almost inconsequential at this stage, as I am reasonably confident in not using such a high-order filter that we hit any timing constraints.

So, I have written some initial matlab code to test my application on a first order Butterworth pair, as seen below. It is worth noting that our sampling frequency is 1 Khz, meaning the nyquist frequency is 500 Hz. As such the normalized cutoff frequency of 0.25 equates to 125 Hz, although this is arbitrary at this stage. Possible problems include the summing of the 2 filters in the last line, or the need to add some integrator term to the high-pass. We also currently have some stock complimentary filters that came with the drone (all came as a big package from Quanser, their autonomous vehicle research suite thingy), but it seems lackluster to just use them and not know what they do. I will also be taking a closer look at the filters for the IMU and magnetometer in a short while after this is complete, as I think I can do a little better. 

clc;
clear;
close all;
w_n=0.25; %cutoff frequency

[zL1,pL1,kL1] = butter(1,w_n,'low');
[numL1,denL1] = zp2tf(zL1,pL1,kL1);
sysL1=tf(numL1,denL1);

[zH1,pH1,kH1] = butter(1,w_n,'high');
[numH1,denH1] = zp2tf(zH1,pH1,kH1);
sysH1=tf(numH1,denH1);

sys1=sysL1+sysH1; 


Making the bode plot of this transfer function gives me something completely different from what I was expecting, as also seen below. 


Any help and guidance is appreciated, many thanks!
You are simply seeing the quantization of the double type. Basically you have a pole and a zero that should cancel each other out. But because of the limited precision of the double data type, it doesn't. So instead of having 0 dB you have a gain of 10^-15 dB give or take.

But in real life,  99.99999% of the time 10^-15 dB = 0 dB

For the record, I think you could use the minreal function in Matlab to cancel the almost equal pole(s) and zero(s). You would then get your unity gain and your flat bode plot with 0 dB and 0 degree phase.

I don't like simply posting links, but a reminder on how floating-point works should help you.

https://en.wikipedia.org/wiki/Floating-point_arithmetic