Freesteel Blog » 2016 » December

Friday, December 30th, 2016 at 9:44 pm - - Flightlogger 1 Comment »

Not been taking time to blog much. I have been lost in a vast pool of non-working projects, and trying to learn scipy and pandas and jupyter, which is an effective replacement for twistcodewiki, which itself was derived from Scraperwiki. The in-browser interactive programming environment Jupyter is everything that Scraperwiki was supposed to have become (eg see data journalism example). But that didn’t happen because every other software engineer I knew assured me in No Uncertain Terms that nobody was ever going suffer the sheer excruciating pain of coding in the browser — when real programmers all want to use vi and git.

But I digress to a world where practically all of my coding is in the browser and I’m using tools like matplotlib at a cycle time of around 10 seconds, which is about a quarter of the time it takes to lose any image file you have generated somewhere on your goddamn disk!

Today I put my new hang-glider data logger (which is held together with lots of glue) into the floor of the elevator and recorded the accelerometer in Z sampling at 100Hz.

azacconly

With my new powers of filtering

b, a = scipy.signal.butter(3, 0.005, 'low')
pZ["fz"] = scipy.signal.lfilter(b, a, pZ.az)

we can get this:

azacconlyf

The art of using sensors is to use a set of them in order to access the underlying reality.

Here is the raw barometer reading as the elevator goes up and down several floors.
baroonly

Once again it’s a good idea to filter it, only this time we set the frequency of that butterworth filter to 0.01 instead of 0.005 because the barometer readings are at 50Hz, and we need the time delay from the filters to be comparable.

baroonlyf

Because it’s nice and smooth after filtering, we can differentiate it by differencing it, using the pF[“df”] = pF.f.diff()*50/10. The times 50 is because it’s at 50Hz, and the divide by 10 is to convert from barometric readings to metres (approximately). The scale is metres per second, and metres per second squared. The elevator travels at about 1.5m/s. The section on the right with the green line lifted for 12 seconds was traveling down about 3 floors at once, so the air pressure was constantly rising.

dbaroacc

Acceleration is the differential of the velocity, so setting ddf = df.diff()*50 gets us the following plot:
ddbaroacc

What we now see are just horrendous oscillations in the barometer. Every time you differentiate you magnify the noise.

But we can reduce the lowpass frequency of this filter by a factor of 5, so it’s 0.001 for the accelerometer and 0.002 for the barometer, and get a curve that agrees (where all the scaling seems to work out without any fudging).

ddbaroaccf

The problem now is you can’t see the elevator moving between the different floors at all. What are those barometric oscillations from? Was I breathing too hard?

Let’s plot the barometer differences against the filtered barometer differences.
dbaronof

There’s a little too much noise, so filter the barometer difference by 0.1 instead and plot it against the 0.01 filter.

dbaronof1

No, I can’t see those oscillations in the underlying data either.

To my eyes they only start to spike up with a filter low pass frequency of 0.03.

dbaronof3

So I’m pretty stumped. But now I’ve got my eye on it, I can see a slight waviness at about this frequency in the original barometric reading in the third picture from the top, so it’s not an illusion or a creation of this filter.

If I hadn’t just glued my entire unit together I might have separated just this device out to see if it’s caused by yet another oscillation in the circuitry, like the kind I’d seen on a temperature sensor.

Maybe I’ll have to rerun this elevator experiment with all the other sensors turned off as much as possible to see if I can get the the bottom of this mysterious 4 second oscillation in the barometer.

Then there’s the temperature decay curves to sort out (from the use of a hairdryer this morning), which are not as good as I had hoped.

I am pleased to have finally implemented an AB Butterworth filter on my arduino.

It’s not a big deal. The code is roughly as follows, and you can get the a[] and b[] arrays from SciPy by executing the function scipy.signal.butter(3, freq, ‘low’):

class ABFilter
{
    float* a, *b, *xybuff; 
    int n, xybuffpos; 
}

ABFilter::ABFilter(const float* la, const float* lb) :
    a(la), b(lb), n(len(la)), xybuffpos(0) // len(a) == len(b)
{  xybuff = (float*)malloc(2*n*sizeof(float));  }

float ABFilter::filt(float x)
{
    xybuff[xybuffpos] = x; 
    int j = xybuffpos; 
    float y = 0; 
    for (int i = 0; i < n; i++) {
        y += xybuff[j]*b[i]; 
        if (i != 0)
            y -= xybuff[j+n]*a[i]; 
        if (j == 0)
            j = n; 
        j -= 1; 
    }
    if (a[0] != 1)
        y /= a[0]; 
    xybuff[xybuffpos+n] = y; 
    xybuffpos += 1;
    if (xybuffpos == n)
        xybuffpos = 0; 
    return y; 
}