Tuesday, June 10, 2014

Weight charts from wwdiary, part 2

In part 1 [1] I described a way of generating nice weight charts from the data stored in the wwdiary [2] weight loss app.

I record a daily weight measurement, usually first thing each morning. As one would expect this can vary a lot from day to day for all sorts of reasons. To get a better handle of the trends a smoothing filter can be applied to the data. There are many methods of doing this, eg Moving Averages. The one I've opted to use here is the Leaky Integrator [3] due to its simplicity.

The unsmoothed chart looks like this.

The following is a GnuPlot chart that applies the Leaky Integrator smoothing to the data (assumed to be in file weight.dat, see part 1 on how to extract this from the wwdiary database):

# GnuPlot script to plot smoothed weight from datafile in format
# uuuuuuuuuu ww.w
# Where uuuuuuuuuu is date in unix epoch time
# and ww.w is weight (in kg here, but trivial to change to other units).
# One record per line. 
# Smoothing is acomplished with the Leaky Integrator
# http://en.wikipedia.org/wiki/Leaky_integrator
# Joe Desbonnet, jdesbonnet@gmail.com, June 2014.
set xdata time
set timefmt "%s"
set format x "%m/%Y"
set xtics font "Arial, 10" 
set title "Smoothed Weight Loss Curve using Leaky Integrator λ=0.8"
set xlabel "Time"
set ylabel "Mass (kg)"
set grid
set key off
set term pngcairo size 800,400
set output "weight.png"
# Lambda (0 ≤ λ < 1.0) determines smoothing. Higher values for more smoothing.
leaky_integrator(x) = (ym1 = ym1==0 ? x : (1-lambda)*x + lambda*ym1)
plot 'weight.dat' using 1:(leaky_integrator($2)) with linespoints, \
     'weight.dat' using 1:2 lt 3

After applying the smoothing filter (and plotting the original points in blue for reference):

You can adjust the 'strength' of the smoothing filter by changing the value of lambda. Lambda can vary from 0 to 1. The closer to 1 the higher the smoothing effect.


[1] http://jdesbonnet.blogspot.ie/2014/06/weight-charts-from-wwdiary-part-1.html
[2] https://play.google.com/store/apps/details?id=com.canofsleep.wwdiary
[3] http://en.wikipedia.org/wiki/Leaky_integrator

Monday, June 9, 2014

Weight charts from wwdiary, part 1

I've been using an excellent Android app called wwdiary [1] for over two years to lose (and then manage) my weight. The app allows you to record your daily weight, but unfortunately has no inbuilt charting [2]. So I wrote a script to extract the data and make my own charts.

Here is my chart:

Here's how to make your own:

This recipe requires a few command line tools. I'm using Linux, but to the best of my knowledge all these tools are available for Windows and Mac and should work exactly the same way described below. You'll need SQLite (there are several implementations), and GnuPlot.

wwdiary uses the SQLite database engine for storage, and you can grab a copy by going into the app settings and saving a copy of the database to SD Card. Copy this file to your computer using whatever file transfer tools you have to hand (I often just email it to myself using the gmail app file-attach option). Let's call this file wwdiary.database.

First you need to extract the weight from the SQLite database in a format that's easy for GnuPlot to parse.  Create a SQLite script file. Let's call it extract_weight.sqlite:

.separator "\t"
select date/1000,weight*0.453592 from weightlog;

Note the "/1000" is to convert the unix epoch timestamp in milliseconds to seconds and the "*0.453592" is to convert imperial pounds used internally in the app [3] to kg. Run the script and send the output to weight.dat:

sqlite3 wwdiary.database < extract_weight.sqlite > weight.dat

The contents of weight.dat should look something like this:

1379890800 76.2999384291992
1379977200 76.0999346446533
1380063600 75.6999409180908
1380150000 75.5999355651855

Create the plot script, call it say weight.gp:

# GnuPlot script to plot weight in format
# uuuuuuuuuu ww.w
# Where uuuuuuuuuu is the unix epoch time in seconds
# and ww.ww is weight (I use kg, but use what you like).
# Joe Desbonnet, jdesbonnet@gmail.com, June 2014.
set xdata time
set timefmt "%s"
set format x "%m/%Y"
set xtics font "Arial, 10" 
set title "Weight Loss Curve"
set xlabel "Time"
set ylabel "Mass (kg)"
set grid
set key off
set term pngcairo size 800,400
set output "weight.png"
plot 'weight.dat' using 1:2 with linespoints

And now run it:

gnuplot weight.gp

All going well, you should find a PNG image file called weight.png in the same folder with your weight chart.

In part 2 I'll show you how to apply a smoothing filter to the chart.


[1] https://play.google.com/store/apps/details?id=com.canofsleep.wwdiary
[2] There is a subscription service available where a chart can be viewed on an online dashboard.
[3] Oh please if you're designing software use SI (metric) units under the hood. Convert to antiquated Olde Worlde/American units at the user interface if you must. Even the Apollo moon avionics software from 1960s USA calculated in SI, converting to feet, furlongs, hundred-weight etc only at the last second at the user interface. (See http://dodlithr.blogspot.ie/2011/10/apollo-guidance-computer-agc-software.html)

Monday, May 26, 2014

Solar FREAKIN' Roadways: snow melting feasibility

Snow removal test (from Solar Roadways
Press Images library).
This video ("Solar FREAKIN' Roadways!") has been circulating on Facebook in the last few days. I find the style rather annoying (the kind that piques my BS-o-meter) and the claim of being able to melt snow off the road hit a major BS alarm bell for me. But rather than poo-poo the idea, I did some back-of-the-envelope calculations on the feasibility of this claim (remember that melting ice takes a lot of energy: the equivalent to heating the same quantity of water by 83.5 degrees C).

I'm going to make the following assumptions in this calculation.
  • Ambient temperature is close to freezing (ie 0C). Ie not that cold, but not warm enough to melt ice.
  • These road tile PV panels are 10% efficient (taking into account losses due to covering material, dirt, and that not all the tile area is covered in PV material, loss due to battery storage, non-optimal angle ie horizontal). This I think this is a very generous assumption.
  • 1m of snow fall melts to 0.1 meters of water [1]
  • Heater is a simple resistive heating element (100% efficient)
So let's assume that 'h' meters of snow needs to be melted over area 'A' (obviously we're dealing with fractions of meters here!)

 (mass  = volume x density of water)

M (kg) = A . h . 1000

The latent heat of fusion of water (334kJ/kg) is the energy required to melt ice from 0C ice to 0C water. So energy requirement is:

E (Joules) = 334x10^6 . A . h

Let's express this energy in terms of sunlight time required. The solar constant (the amount of energy arriving per  unit area that would be incident on a plane perpendicular to the rays) is 1kW/m^2. So with 10% efficiency that's 100W / m^2

So sunlight time needed is t = E / 100W

t (seconds) = 334x10^4 . A . h

So for 1cm of snow, translating to 1x10^-3 m of water and expressing on a per m^2 (ie A=1)
t  = 3340 seconds. A little under one hour.

So subject to the very generous assumptions above, one hour of stored sunshine can melt 1cm of snow.

Now bear in mind that snowfall occurs during overcast periods, so we're going to have to rely on storage of this sunlight energy (presumably in electrochemical batteries, although there may be other ways, eg phase change materials). Looking at the current (2014) costs of Lithium ion rechargeable batteries, I'm seeing figures around $500 (USD) per kWh of storage [2][3][4]. So $500 in battery costs gives you about 1m of snow melting capability on a 1 meter square panel.

So the idea isn't complete bunkham. It's within the realm of physical possibility. This will only work in certain places where snowfall is relatively light and intermixed with regular spells of sunshine. Obviously in areas where the sun does not shine (canyons, forests, extreme latitudes) it just won't work. It takes just one heavy snowfall to exhaust the batteries and after that they'll be starved of sunlight until the snow melts naturally. Also remember that snow falls in winter when the sun is low in the sky and only for a few hours a day.


I notice the heating aspects of these tiles is briefly discussed (but without any hard numbers) on their website FAQ here [5] and general information about Solar Roadways is here [6].

Dave Jones of the EEVBlog has a rather damning analysis of the power budget for the LED lighting (let alone the heating aspect) here [7]


[1] Snow-water equivalent
[2] http://www.forbes.com/sites/sarwantsingh/2014/02/27/elon-musks-500000-unit-battery-plant-myth-or-reality/
[3] http://www.mckinsey.com/insights/energy_resources_materials/battery_technology_charges_ahead
[4] http://en.wikipedia.org/wiki/Electric_vehicle_battery#Battery_cost_and_parity
[5] http://www.solarroadways.com/faq.shtml#faqHeat
[6] http://www.solarroadways.com/intro.shtml
[7] https://www.youtube.com/watch?v=obS6TUVSZds


26 May 2014: My initial figure of $160 (USD) per kWh of battery storage (2014 prices) seems to be very much on the low side. Seeing more consensus at $500 per kWh. See references [2][3].
26 May 2014: Added links to the company website and a link to a FAQ where the heating of the tiles is discussed. Also image of snow remove test added (taken from Press Images page.
24 June 2014: Added link to Dave Jone's analysis of the project

Note to commenters: I've got moderation enabled, so please allow a day or two for your comment to appear here.

Tuesday, February 18, 2014

Capturing a thermal IR video with ffmpeg on Linux

When tethered to a computer with USB the Flir E4 [1] thermal camera can act as webcam. The frame rate is limited to 9fps due to ITAR restrictions [2], but this is sufficient for most applications. On Linux FFMPEG [3] can be used to capture a video with this command line:

ffmpeg -f v4l2 -r 9 -s 640x480 -i /dev/video0 thermal.avi

You can preview while shooting with the camera's LCD. But sometimes you want to see a preview on screen. That turns out to be tricky.

I found a useful solution on this thread [4]:

1. Create a fifo file:
mkfifo /tmp/livevideo.fifo

2. Run mplayer on the fifo:
mplayer  -cache 512 -really-quiet /tmp/livevideo.fifo

3. Run ffmpeg to capture to file and also relay to the fifo:
ffmpeg -f v4l2 -r 9 -s 640x480 -i /dev/video0 thermal.avi -map 0 -f avi -vcodec rawvideo -y /tmp/livevideo.fifo

Here is a sample video captured with ffmpeg:


[1] also the Flir E6, E8, the Exx range and most (all?) recent models as far as I know

[2] http://en.wikipedia.org/wiki/International_Traffic_in_Arms_Regulations

[3] http://www.ffmpeg.org/

[4] http://nerdfever.com/?p=2661

Monday, February 17, 2014

Sox spectrogram log frequency axis and upper/lower frequency limits

Linear spectrogram
As a result of some of the work I did last week on rendering MP3 files to scrolling spectrum waterfall plots, I noticed that most of the interesting detail was in frequencies under 1kHz.  One obvious way to see the detail at lower frequencies, yet still keep the higher frequencies is to plot on a log frequency scale.

Unfortunately this is currently not supported in SoX, so I modified the spectrogram module to add this feature. This I noticed that plotting charts from 1Hz to the nyquist frequency meant that a lot of screen space was wasted in the lower 1 - 50Hz range which is not the most interesting for most audio files. So I also added a switch to set the lower and upper frequencies of the chart (this actually took considerably more work than the plotting to the log scale!).

Log axis spectrogram. Lower frequency details are far more visible.
The modified spectrogram.c file is available here [1]. Visit the SoX site [2], download the source code, replace src/spectrogram.c with this file and compile. This has been tested with the latest code as of 14 Feb 2014.

Summary of changes:
  • Implement -L which plots the spectrogram on a log10 frequency axis. By default from 1Hz to nyquist frequency.
  • Implement -R <low_freq>:<high_freq> : restrict spectrogram frequencies. I also updated the linear scale code to honor this switch.


sox mymusic.mp3 -n spectrogram -L -R 50:8k

Known problems and observations:
  • At lower frequencies (1 - 100Hz) each spectrum bin in q->dBfs[] array occupies several pixel rows: so it looks blocky. A lesser problem: at the top some frequency bins will be ignored because there will be more than one bin for each pixel row. I could fix the blockyness by getting a more detailed spectrum.
  • If lower and upper frequencies are in the same decade then there are no y axis tick labels. Hopefully will have a fix for that soon.
  • When compiling I'm getting this warning when using log10f() and powf() : "warning: passing argument 1 of 'log10f' as 'float' rather than 'double' due to prototype". According to the man pages for those function they should accept float args! I see some discussion on this relating to compiler switches: but I don't want to go changing anything there.
  • I created two new functions to parse the -R frequency range switch to cut down on code duplication.
    parse_range (const char *s, int *a, int *b)
    parse_num_with_suffix (const char *s, int *a)
    I'm not familiar with the sox code base, so I'm not sure if similar functions already exist (I couldn't find any). I'm also concerned I'm making a temporary change to a string in parse_range when the type is const char *.  I need to brush up on my C :)


[1] https://github.com/jdesbonnet/joe-desbonnet-blog/tree/master/projects/sox-log-spectrogram

[2] http://sox.sourceforge.net/

Saturday, February 15, 2014

Boosting audio volume in a video file

For some reason I'm finding that some video files played directly on my new Samsung TV (by mounting a USB drive) have a low volume and boosting the TV's volume beyond 40 (of presumably 100) yields no additional gain.

This ffmpeg command can be used to boost the volume in the file:

ffmpeg -i myvideo.mkv -vcodec copy -af "volume=12dB" myvideo-boosted.mkv

You can vary the boost by chaning the volume gain parameter. I found 12dB was about right. Negative values are allowed also.

Wednesday, February 12, 2014

Convert MP3 to a scrolling spectrum waterfall plot video

There are many utilities that display a scrolling spectrum waterfall plot [1] from a signal, but I was unable to find any open source utility that converted an audio file into a video file with a scrolling waterfall plot + the audio.

The SoX [2] sound utility can generate a static spectrum waterfall plot image from an audio file (or part of it), but it can't make a video. So I wrote a script to do this.

It's a very brute force approach. It requires lots of CPU time and lots of temporary disk space. The script depends on SoX, GNU Parallel, mencoder (or ffmpeg). GNU Parallel is optional, but will result in significant speed up on a multi-core system.

To use, create an empty directory on a volume with plenty of disk space and run with the audio file as a parameter.  Eg:

./make-spectrogram-video.sh -t "My Music File" mymusic.mp3

The output will be written to output.avi and output.mp4. Other options include setting frame rate, the speed of the scrolling, audio credit text etc. To get full help do this:

./make-spectrogram-video.sh -h

The script is available on GitHub here [2]. Here is a sample output video of Bach's Toccata and Fugue in D Minor [3] :


17 Feb 2014: I noticed that all the interesting details in music is squashed down at the very bottom of the spectrogram. So I updated the spectrogram module in SoX to have the option of plotting on a log axis. This isn't in the offical sox distribution. See this blog post for more details [4].


[1] http://en.wikipedia.org/wiki/Waterfall_plot

[2] https://github.com/jdesbonnet/audio-to-waterfall-plot-video

[3] Music MP3 file from https://archive.org/details/ToccataAndFugueInDMinor. YouTube video at http://www.youtube.com/watch?v=utp95bprqeg

[4] http://jdesbonnet.blogspot.ie/2014/02/sox-spectrogram-log-frequency-axis-and.html

Tuesday, February 11, 2014

Simple is better when it comes to microwave oven UI

About a year ago the magnetron (the magic yolk that generates microwave radiation) on my old Sharp microwave oven failed. This is not economically repairable failure mode, so I went shopping for a new oven. I resisted the temptation to get one with lots of lights, buttons and bluetooth. Instead I opted for a €49 oven with a simple mechanical user interface.

After a year of use I have to say it's the best oven user interface I've ever used. It's faster to set than the keypad of the old oven. The only feature I miss from the old oven is a kitchen timer. And it's a missed opportunity here: all they needed to do was add a zero power option on the power setting knob!

Monday, February 3, 2014

A fine stretch in the evenings

This is small web app (click to try) I wrote, inspired by that phrase often said at this time of the year at this latitude: "a fine stretch in the evenings". It is also inspired by a low-latitude friend who thought our long summer days a bit freaky :-)  [1]

It's a well known phenomena that day/night lengths vary throughout the year with latitude. From no variation at the equator to the extreme 24 hour days and nights near the poles. This app helps you visualize this variation. Enter your latitude/longitude or choose your location on a map. You'll get a visualization of the length of days/nights throughout the year at that location. The system will estimate the appropriate timezone for each location and adjust for 'wall clock' time accordingly.

On the chart black is night (obviously), yellow is day and the blue band is twilight. If there is daylight saving in the timezone then you will see a discontinuity in the sun rise/set curve during the transitions.

You can have one or two locations (although a HD monitor will be required to fit both charts side by side).

Some interesting well-known observations:

1. Here in Ireland the days get depressingly short during the winter and we all get a bit cranky, sick, moldy and drink too much. [2]

2. At the equator, day/night is pretty much same all year around. They all have a happy disposition and they don't drink at all! [3]

3. Up near the arctic circle there comes a point in the year where the sun doesn't rise at all. Oddly up there there they don't seem to get nearly as cranky, sick or moldy as we do in Ireland, but they still drink too much.

Some lesser known observations:

1. The length of twilight varies with latitude. In fact you can estimate latitude using just that [4] 

2. The shortest day of the year, the day of latest sun rise and the day of earliest sun set are not the same! (ie the sun rise and set curves are not exact mirrors of each other)

This app uses HTML5 has been tested with the latest Chrome and Firefox browsers. It will probably will cause Internet Explorer to explode.
[3] There are exceptions!

[4]  http://www.birdtracker.co.uk/ 

Friday, January 31, 2014

Generating a spectrogram of The Shepard Tone

Earlier today I saw a tweet from Tim O'Reilly pointing to an interesting article in The Atlantic on audio illusions.  The Shepard Tone illusion sounds like a continuously increasing or decreasing tone which goes on indefinitely... which is impossible because you'd eventually end up at at inaudibly low or high pitch.

Curious as to what was going on I downloaded the YouTube video, extracted the audio and ran off a spectrogram with the open source sox audio tool. I had done this before last year, so I was familiar with the process.

1. Download the YouTube video:
youtube-dl -o shepard_tone.mp4 http://www.youtube.com/watch?v=DfJa3IC1txI

2. Extract the audio from the video
ffmpeg -i shepard_tone.mp4 -f mp3 shepard_tone.mp3

3. Convert to mono (the YouTube video has stereo audio)
sox shepard_tone.mp3 shepard_tone_mono.mp3 remix 1,2

4. Create the spectrogram
sox shepard_tone_mono.mp3 -n spectrogram -o spectrum.png

One small problem: most of the action is down at the lower frequencies. So I resampled the audio and regenerated the spectrogram so that the interesting bits were more visible.

5. sox shepard_tone_mono.mp3 -r 4k -o shepard_tone_4k.wav

6. sox shepard_tone_4k.wav -n spectrogram -o spectrum.png

Tools used :