## Revisiting Hum Filters

I’m in a relatively noisy environment when it comes to mains hum having overhead power lines nearby. So any ELF/VLF receiver I put together will have to deal with this.

To get an idea of what kind of filtering I might need I simply put a jack plug with bare terminals into the mic in of a laptop, held on to it to make myself an antenna, and recorded for half a minute using Audacity. Here’s a snippet of the resulting waveform:

Yup, that is one well-distorted sine wave. Reminds me of the waveform going to bulbs from triac-based dimmers, though haven’t any in the house.

More usefully, here’s the spectrum plot (Audacity rocks!) :

There’s a clear peak at 50Hz. Next highest is at 150Hz, the 3rd harmonic. It’s around 12dB down, which (assuming it’s the voltage ratio being shown, ie. 20*log10(V2/V1)) is 1/4 of the voltage. Next comes 100Hz, the 2nd harmonic, about 30dB down, about 1/32 of the voltage (from ratio = 10^(dB/20)).

(I’m in Italy where like most of the world the mains AC frequency is 50Hz. In the Americas it tends to be 60Hz).

So I reckon I definitely need to cut the 50Hz as much as possible, probably 150Hz too.

Digital filters are relatively straightforward to implement in software, but here there’s a snag. The incoming signal is analog, so will need to go through an ADC. The ELF/VLF signal of interest is likely to be of very small amplitude compared to the mains hum. So using say a 16-bit ADC, capturing the whole signal at maximum resolution, it’s conceivable that the interesting signal only occupies a couple of bits, or maybe even be below a single bit. So really the filtering has to happen in the analog chain, before the ADC. Experimentation will be needed, but I imagine a setup like this will be required:

There are a few options for the kind of receiver to use, essentially coil-based (magnetic component of the radio wave) or antenna (electrical component), the nature of the early circuitry and pre-amp will be dependent on this. But the main role of the pre-amp is to boost the signal well above the noise floor of subsequent stages (using low-noise components). At a first guess, something in the region of x10 – x100 should be adequate.

Next comes the filter(s). Now the fortunate thing here is that the ELF/VLF frequency ranges I’m considering, say 5Hz-20kHz are pretty much the audio frequency ranges and are thus within the scope of standard audio components. Well, 5Hz is below the nominal 20Hz-20kHz figures given for audio, but the key thing is that at the high end, it’s nowhere near anything requiring exotic components. Even the humble 741 op-amp (dating from 1968) has a unity-gain bandwidth around 1MHz. For the TL071 family, a reasonable low-cost default these days it’s 3MHz.

One option for filtering the mains hum out would be to use a high pass filter and only look at the higher end of VLF (conversely, a low pass filter and go for ELF). But notch (band stop) filters can be pretty straightforward, so it should be productive to target just the 50Hz (and maybe 150Hz).

(A more exotic approach would be to use something like an analog bucket brigade line device as used in many analog phaser & flanger audio effects boxes, with its delay fixed at the period of the fundamental 50Hz. Mixing this inverted with the input signal non-inverted will cause cancellation at the fundamental and all it’s harmonics, ie. a comb filter. But not only does that seem overkill here, it will in effect degrade the signal of interest).

There are a few alternatives for notch filters. While they can be built from passive components, there are significant benefits to using active components, especially in terms of controlling the parameters. For these reasons and circuit simplicity, op-amps are a good choice over discrete components.

There are three leading candidate circuit topologies, as follows.

#### Active Twin-T Notch

This classic passive circuit is the starting point.

The notch frequency is given by fc = 1 / (2 pi R C)

This assumes a low impedance source for Vin and a high impedance connected to Vo, which can easily be achieved using op-amp buffers. One drawback of this setup is that its selectivity, the slope of the sides of the notch, is fairly poor. This can be significantly increased by using op-amps to bootstrap the T :

The notch frequency is determined as for the grounded T above, only this time the Q/selectivity can be varied, according to the values of R4 and R5.

But a troublesome problem remains: all 6 components on which the frequency depends have to have precise values to place the notch where required. Any variation is likely to lead to a sloppy notch, of low Q. While 1% tolerance resistors are the norm these days, capacitors tend to have tolerances more like 5 or 10%.  One option is to use reasonably well-matched capacitors (from the same batch) and vary the resistors. But this still leave 3 variables, with some level of interdependence.

(I’ve actually got this one on a breadboard at the moment. For a one-off circuit it isn’t unreasonable to use resistors a little below the calculated values in series with pots, and once fine tuned replaced with fixed values).

#### Bainter Notch

This is quite a nifty circuit (and new to me). The main benefits are described in the title of Bainter’s own description : Active filter has stable notch, and response can be regulated. Notch depth depends on gain and not (passive) component values.

I’ve yet to play with this one, but it certainly shows promise. A downside is that the component values calculation is rather unwieldy. Another ref. is this TI doc: Bandstop filters and the Bainter topology.

#### State Variable Filter

The State Variable topology is very versatile, offering high- and low-pass outputs as well as bandpass. By mixing the high- and low-pass outputs or the input with the bandpass output, a notch can be achieved. Crucially the gain, center frequency, and Q may be adjusted separately. A bonus compared to the Twin-T is that the frequency is determined by just 2 resistors and 2 capacitors. A few days ago I stumbled on a tweaked version of the standard topology which offers a few advantages. I won’t go into details here, it’s all described in the source of this diagram – Three-op-amp state-variable filter perfects the notch.

Once I’ve played with the Twin-T a bit more, I’ll have a go with this one. I have a good feeling about it.

## Hall Effect-based Seismometer, Sanity Check Experiments

PS. Oops! I made a silly mistake in the breadboarding, if you look closely at the photo you can see that the 10k ground resistor at the input of the op amp is going to + input, not – as intended. Which kind of messes up all my measurements. Hey ho. I have since made a ball-bearing in a jar (1 axis) sensor and roughed out a signal conditioning circuit (which will now need tweaking…), so will repeat the experiment here and do another post asap.

A fun part of this project is the investigation of hardware possibilities for detecting seismic events and ELF/VLF signals. Even though I’m aiming towards minimum budget hardware, my funds for this have been virtually non-existent so I’ve not got much done (grumble, grumble).

For a seismometer, the requirements as it seems to me, are: simplicity, reasonable sensitivity and low cost. Ideally I want to monitor all 3 dimensions with relatively wide bandwidth. A non-requirement is any kind of absolute accuracy or calibrated measurement.

There are a variety of options for seismic sensors, most that I’ve seen fall down for these requirements in one way or another. I won’t go into them here – try searching for accelerometers (low sensitivity), geophones (expensive), pendulum-based systems (complicated build, 3 dimensions would be very tricky…). To give a ballpark, prices for a ready-made seismometer system based on the Raspberry Pi, the Raspberry Shake, start at \$375 USD. That’s for one dimension, using a geophone sensor.

Almost a year ago I sketched out an idea for something that might work.

At the time I picked up a linear Hall Effect device from Jaycar, a UGN3503UA, costing just \$7.75 AUS. It’s in case very like a transistor, just 3 pins : +ve, -ve power and output. An example use in the datasheet uses the same principle as I want to exploit:

A magnet is glued to the back of the sensor. As a (ferrous material) cog approaches the sensor, the magnetic field increases, correspondingly increasing the devices output voltage.

The other day a bag of ball bearings arrived. I just got around to having a play. This is what the setup looked like:

I’ve got the Hall Effect device soldered to a connector to make breadboarding easier. On the left of it is a blob of Blue Tack attaching a 1cm diameter/3mm deep neodymium magnet. On the right, a 5/8″ steel ball precision mounted between my finger & thumb.

Right now I’ve only got a crude +/- homebrew power supply, so I’m using an op amp to buffer a potential divider to provide a lower voltage to suite the device. Another op amp is used to provide a 10x amplifier from the output of the device.

When I put the magnet in direct contact with the sensor it saturated it at one extreme or the other. I seemed to get best results with around 1cm space in between. With a 5.2v supply to the sensor, this led to a no-magnet output of 2.52v (after the 10x amplification). With the magnet, this changed to 3.07v or 1.76v depending on polarity. With the ball bearing at 1cm away this changed by approx 0.01/0.02v, steadily increasing from there to 3.50/1.22v when the ball bearing touched the sensor.

This sensitivity was less than I’d hoped, but will hopefully be enough to be usable if I tweak a few of the components. I reckon it’s definitely worth going for a prototype, see how it behaves in practice.

I’ll need to find a very small jar 🙂

Here are my full notes:

## Human Impact on Radio Nature

I’ve stumbled on two pieces of info related to this in the past couple of days so reckon it’s worth making a note.

The first is NASA’s Van Allen Probes Spot Man-Made Barrier Shrouding Earth, actually about high-powered VLF transmitters for ground-submarine comms,  probably affecting the near-space environment. “A number of experiments and observations have figured out that, under the right conditions, radio communications signals in the VLF frequency range can in fact affect the properties of the high-energy radiation environment around the Earth”. The main reference paper is on a pay-for site, so little detail is at hand.

The second is Why I Quit Natural Radio, a post by a Radio Nature enthusiast who’s been monitoring ELV/VLF for decades, noting a massive drop off in the ‘interesting’ natural signals he receives. He suggests the cause may be the rise in the use of mobile phones, associated UHF/microwave emissions (‘Cellular frequencies‘) affecting the magnetosphere and/or ionosphere thus impacting VLF propagation. A term he’s coined is rather disturbing ‘electromagnetic smog’.

He also refers to HAARP, a research system (and favourite of conspiracy theorists) that has historically blasted the ionosphere with high power HF. According to official sources it hasn’t been used for a long time. I have my doubts that relatively brief, localised high-energy signals like these would have any lasting impact – similar events might well occur in nature, and these natural systems tend to be very resilient.

## Preconditioning Seismic Data

The filtered data I have is CSV with lots of lines with the fields :

`datetime, latitude, longitude, depth, magnitude`

The latter 4 fields will slot in as they are, but a characteristic of seismic events is that they can occur at any time. Say today 4 events were detected at the following times:

E1 01:15:07 lat1 long1 d1 2.2
E2 01:18:06 lat2 long2 d2 3.1
E3 01:20:05 lat3 long3 d3 2.1
E4 08:15:04 lat4 long4 d4 3.5

To get the data in a shape that can act as input to a neural network (my first candidate is PredNet), it seems like there are two main options:

#### Time Windows

Say we decide on a 6 hour window starting at 00:00. Then E1, E2, E3 will fall in one window, E4 the next.  Which leads to the question of how to aggregate the first 3 events. Often events are geographically clustered, a large event will be associated with nearby foreshocks and aftershocks. For a first stab at this, it doesn’t seem unreasonable to assume such clustering will be the typical case. With this assumption, the data collapses down to :

[00:00-06:00] E2 lat2 long2 d2 3.1
[06:00-12:00] E4 lat4 long4 d4 3.5

This is lossy, so if say E1 and E2 were in totally different locations, the potentially useful information of E1 would be lost. A more sophisticated strategy would be to look for local clustering – not difficult in itself (check Euclidian distances), but then the question would be how to squeeze several event clusters into one time slot. As it stands it’s a simple strategy, and worth a try I reckon.

#### Time Differences

This strategy would involve a little transformation, like so:

E1[datetime]-E0[datetime] = ? lat1 long1 d1 2.2
E2[datetime]-E1[datetime] = 00:03:01 lat2 long2 d2 3.1
E3[datetime]-E2[datetime] = 00:02:01 lat3 long3 d3 2.1
E4[datetime]-E3[datetime] = 07:05:01 lat4 long4 d4 3.5

Now I must confess I really don’t know how much sense this makes, but it is capturing all the information, so it might just work. Again, it’s pretty simple and also worth a try.

I’d very much welcome comments and suggestions on this – do these strategies make sense? Are there any others that might be worth a try?

## Seismic Data – fixed?

As described in my last post, I was seeing significant gaps in the seismic event data I was retrieving from the INGV service. So I re-read their docs. Silly me, I’d missed the option to include query arguments restricting the geo area of the events (I had code in a subsequent script doing this).

While tweaking the code to cover these parameters I also spotted a really clumsy mistake. I had a function doing more or less this –

```for each event element in XML DOM:
extract event data
return list```

D’oh! Should have been –

```for each event element in XML DOM:
extract event data
return list```

I’ve also improved error handling considerably, discriminating between genuine HTTP errors and HTTP 204 No Content. Now I’ve narrowed the geo area and reduced the time window for each GET down to 1 hour, there are quite a lot of 204s.

I’m now running it over the time period around the l’Aquila quakes as a sanity check. Jeez, 20+ events in some hours, 10+ in most.

Assuming this works ok, I’ll run it over the whole 1997-2017 preiod, hopefully in ~12 hours time I’ll have some usable data.

PS. Looking good, for the 30 days following that of the l’Aquila big one, it produced:

in_zone_count = 8877
max_depth = 62800.0
max_magnitude = 6.1

## Seismic Data Wrangling

Following my interim plan of proceeding software-only (until I’ve the funds to get back to playing with hardware), I’ve been looking at getting seismic event data from the INGV Web Service into a Keras/Tensorflow implementation of PredNet.

My code is on GitHub, and rather than linking to individual files which I may rename, I’ll put a README over there with pointers.

As a first step, I put together code to pull of the data and dump it down to simple CSV files. This appeared to be working. The demo implementation of PredNet takes HDF5 data from the KITTI  vision dataset (videos from a car on road around Karlsruhe), extracting it into numpy arrays, with the PredNet engine using Keras. To keep things simple I wanted to follow the same approach. I’m totally new to HDF5 so pinged Bill Lotter of the PredNet project for clues. He kindly gave me some helpful tips, and concurred with what I’d been thinking – keep the CSV data, process that into something PredNet can consume.

The data offered by the Web Service is good XML delivered over HTTP (props to INGV). But it does include a lot of material (provenance, measurement accuracy etc) that isn’t needed here. So my service-to-CSV code parses out just the relevant parts, producing a line for each event:

```datetime, latitude, longitude, depth, magnitude

e.g.

2007-01-02T05:28:38.870000, 43.612, 12.493, 7700, 1.7```

I couldn’t find the info anywhere, but it appears that the INGV service records go back at least to somewhere in the early 1990’s, so I chose 1997-01-01T00:00:00 as a convenient start datetime, giving me 20 years of events.

For this to be a similar shape to what PredNet expects, I will aggregate events within a particular time period (actually I think taking the most significant event in that period). I reckon 6 hour periods should be about right. This also seemed a reasonable window for calling the service (not). I’ll filter down the events to just those within the region of interest (northern Italy, see earlier post)  and then scale the latitude & longitude to an easy integer range (probably [128, 128]). For a first pass I’ll ignore the depth field.

As it happens, I’m well on the way to having implemented this. But along the way I did a few sanity checks, eg. checking for maximum event magnitude in the region of interest, (I got 4.1), and it turned out I was missing some rather significant data points.  Here’s one I checked for:

The 2009 L’Aquila earthquake occurred in the region of Abruzzo, in central Italy. The main shock occurred at 03:32 CEST (01:32 UTC) on 6 April 2009, and was rated 5.8 or 5.9 on the Richter magnitude scale and 6.3 on the moment magnitude scale; its epicentre was near L’Aquila, the capital of Abruzzo, which together with surrounding villages suffered most damage.

Nope, it wasn’t in the CSV, but the Web Service knows all about it:

http://webservices.ingv.it/fdsnws/event/1/query?eventId=1895389

Doing a service call for that whole day:

http://webservices.ingv.it/fdsnws/event/1/query?starttime=2009-04-06T00:00:00&endtime=2009-04-06T23:59:59

–  yields 877 events – nightmare day!

I’d set the timeout on the HTTP calls to 2 seconds, but there is so much data associated with each event that this was woefully inadequate. Since upped to 5 mins.

Manually checking calls, I was also sometimes getting a HTTP Status Code of 413 Request Entity Too Large. This puzzled me mightily – still does actually. It says request entity, not requested (or response) entity, but the way it’s behaving is that the response requested is too large. Either way I reckon the spec (latest is RFC7231) is a little open to misinterpretation here. (What the heck – I’ve mailed the IEFT HTTP list about it – heh, well well, I’ve co-chaired something with the chair…).

Anyhow, I’ve also tweaked the code to make calls over just 1 hour windows, hopefully it’ll now get the stuff it was missing.

Hmm…I’ve got it running now and it’s giving errors throughout the year 2000, which should be trouble-free. I think I’ll have to have it make several passes/retries to ensure I get the maximum data available.

Drat! It’s giving me Entity Too Large with just 1 hour windows, e.g.

http://webservices.ingv.it/fdsnws/event/1/query?starttime=2000-12-13T01:00:00&endtime=2000-12-13T02:00:00

I need to fix this…

## Candidate Neural Network Architecture : PredNet

While I sketched out a provisional idea of how I reckoned the network could look, I’m doing what I can to avoid reinventing the wheel. As it happens there’s a Deep Learning problem with implemented solutions that I believe is close enough to the earthquake prediction problem to make a good starting point : predicting the next frame(s) in a video. You train the network on a load of sample video data, then at runtime give it a short sequence and let it figure out what happens next.

This may seem a bit random, but I think I have good justification. The kind of videos people have been working with are things like human movement or motion of a car. (Well, I’ve seen one notable, fun, exception : Adversarial Video Generation is applied to the activities of Mrs. Pac-Man). In other words, a projection of objects obeying what is essentially Newtonian physics. Presumably seismic events follow the same kind of model. As mention in my last post, I’m currently planning on using online data that places seismic events on a map – providing the following: event time, latitude, longitude, depth and magnitude. The video prediction nets generally operate over time on x, y with R, G, B for colour. Quite a similar shape of data.

So I had a little trawl of what was out there.  There are a surprisingly wide variety of strategies, but one in particular caught my eye : PredNet. This is described in the paper Deep Predictive Coding Networks for Video Prediction and Unsupervised Learning (William Lotter, Gabriel Kreiman & David Cox from Harvard) and has supporting code etc. on GitHub. Several things about it appealed to me. It’s quite an elegant conceptual structure, which translates in practice into a mix of convnets/RNNs, not too far from what I anticipated needing for this application. This (from the paper) might give you an idea:

Another plus from my point of view was that the demo code is written using Keras on Tensorflow, exactly what I was intending to use.

Yesterday I had a go at getting it running.  Right away I hit a snag: I’ve got this laptop set up for Tensorflow etc. on Python 3, but the code uses hickle.py, which uses Python 2. I didn’t want to risk messing up my current setup (took ages to get working) so had a go at setting up a Docker container – Tensorflow has an image. Day-long story short, something wasn’t quite right. I suspect the issues I had related to nvidia-docker, needed to run on GPUs.

Earlier today I decided to have a look at what would be needed to get the PredNet code Python3-friendly. Running kitti-train.py (Kitti is the demo data set) led straight to an error in hickle.py. Nothing to lose, had a look. “Hickle is a HDF5 based clone of Pickle, with a twist. Instead of serializing to a pickle file, Hickle dumps to a HDF5 file.“. There is a note saying there’s Python3 support in progress, but the cause of the error turned out to be –

`if isinstance(f, file):`

file isn’t a thing in Python3. But kitti-train.py was only passing a filename to this, via data-utils.py, so I just commented out the lines associated with the isinstance. (I guess I should fix it properly, feed back to Hickle’s developer.)

It worked! Well, at least for kitti-train.py. I’ve got it running in the background as I type. This laptop only has a very wimpy GPU (GeForce 920M) and it took a couple of tweaks to prevent near-immediate out of memory errors:

```export TF_CUDNN_WORKSPACE_LIMIT_IN_MB=100

kitty_train.py, line 35
batch_size = 2 #was 4```

It’s taken about an hour to get to epoch 2/150, but I did renice Python way down so I could get on with other things.

#### Seismic Data

I’ve also spent a couple of hours on the (seismic) data-collecting code. I’d foolishly started coding around this using Javascript/node, simply because it was the last language I’d done anything similar with. I’ve got very close to having it gather & filter blocks of from the INGV service and dump to (csv) file. But I reckon I’ll just ditch that and recode it in Python, so I can dump to HDF5 directly – it does seem a popular format around the Deep Learning community.

Yes, that to think about too.

My gut feeling is that applying Deep Learning to the seismic data alone is likely to be somewhat useful for predictions. From what I’ve read, the current approaches being taken (in Italy at least) are effectively along these lines, leaning towards traditional statistical techniques. No doubt some folks are applying Deep Learning to the problem. But I’m hoping that bringing in radio precursors will make a major difference in prediction accuracy.

So far I have in mind generating spectrograms from the VLF/ELF signals. Which gives a series of images…sound familiar? However, I suspect that there won’t be quantitatively all that much information coming from this source (though qualitatively, I’m assuming vital).  As a provisional plan I’m thinking of pushing it through a few convnet/pooling layers to get the dimensionality way down, then adding that data as another input to the  PredNet.

Epoch 3/150 – woo-hoo!

#### PS.

It was taking way too long for my patience, so I changed the parameters a bit more:

```nb_epoch = 50 # was 150
batch_size = 2 # was 4
samples_per_epoch = 250 # was 500
N_seq_val = 100 # number of sequences to use for validation```

It took ~20 hours to train. For kitti_evaluate.py, it has produced some results, but also exited with an error code. Am a bit too tired to look into it now, but am very pleased to get a bunch of these: