Sunday, April 13, 2014

Automatic detection of MH370 ULB pings

In my earlier post (Analysis of suspected MH370 pings) I manually inspected a audio waveform to determine ping times of arrival. Whilst this is managable for the short segment I looked at, to perform this analysis over many hours or days of recordings would be a bit tedious and probably error prone.

This afternoon I put together a quick MATLAB script to automate this process and ran it over the whole Youtube video audio.

The algorithm is as follows:
  • Read acoustic data from file in 1s segments
  • For each 1s segment use a Hilbert Transform (using FFT) to calculate the envelope of the waveform
  • Calculate noise by smoothing the envelope with a window of length 0.1s (or 10x expected pulse width)
Envelope (red); Smoothed (black)
  • Calculate SNR as the envelope amplitude divided by the noise
  • Locate rising edges in the SNR that exceed the detection threshold (ie start of pulse)
SNR (green) with detected pulse (x)
  • At end of processing plot detection times of arrival module ping interval (1.106s from previous analysis)
Time of Arrival Ping Graph


In the resulting plot that there are 2 distinct sections of coherent detections. The others are false detects from other in-band transients.

Interestingly I had not noticed the 2 individual segments previously, but upon a second listen a distinct discontinuity can be heard hence the 'realignment' of the ping times of arrival. This plot shows that the first segment corresponds to a time where the TPL was going away from the source as the times of arrival are increasing. The second segment has the TPL coming closer to the source near CPA due to the decrease and stabilisation of time of arrival.

Wednesday, April 9, 2014

Analysis of suspected MH370 pings recorded by Ocean Shield

So this is a bit of a deviation from my normal MATLAB oriented posts. In-fact MATLAB was not used once in this analysis!

I've been following the MH370 disappearance since it happened. With such a lack of confirmed information there have been a lot of interesting theories discussed on multiple forums (Whirlpool and PPRUNE are a couple of interesting ones).

As I have a knowledge in acoustics (I am an Analyst / Software Engineer for an underwater acoustics company and I specialise in tracking things with pingers fitted to them), I've been trying to spread facts when I know them and try and dispel myths.

Today I saw an interesting video from Sky News on a forum of the suspected MH370 pings recorded by ADV Ocean Shield using the SUBSALV Towed Pinger Locator TPL-25 (pinger audible 18 seconds in):


I thought I would see if I could glean any information from the audio in the publicly released video. In particular I was wanting to look at the inter-ping interval to see how steady it was to try and determine if the noise is man-made or natural (ie marine mammal noise). 
The signal is quite noisy with ambient ship noise and other transients (scuffing of shoes on deck!) but a repeated click can be heard.

Process

  • Extracted the audio from the Youtube video using a web-based MP3 extractor
  • Opened the MP3 in Adobe Audition and trimmed the file to the time segment corresponding to good audio (start time 41.787s from beginning of file, duration 22.462 seconds)
  • Bandpass filtered the audio between 2kHz and 10kHz to remove LF ambient noise (NOTE: the 'pinger' is audible at about 3.4kHz however it obviously has been demodulated down from the raw signal observed on the TPL-25 to bring into the audible range)
Filtered Time Domain Signal - Red ticks added at repeated 'ping' noise
  • A manual search of the waveform (zooming in tight where necessary) was performed to measure the start time of each ping transient to within a couple of wavelengths

Results

The detection times of the pings are listed in the table below along with some basic statistics.

Ping #
Time (s)
Interval (s)
1
0.758
2
1.861
1.103
3
2.964
1.103
4
4.067
1.103
5
5.171
1.104
6
6.274
1.103
7
7.38
1.106
8
8.486
1.106
9
9.592
1.106
10
10.696
1.104
11
11.801
1.105
12
12.907
1.106
13
14.014
1.107
14
15.12
1.106
15
16.226
1.106
16
17.333
1.107
18
19.549
1.108
19
20.657
1.108
20
21.764
1.107
Mean
1.10544
Std Dev
0.00172

Single Ping zoom view - Selected portion is 12ms duration

 As you can see, the ping interval is very stable with an average of 1.105s and standard deviation of less than 2ms.

This is incredibly stable over the 20 seconds that I've looked at and very unlikely to be a marine mammal (ie beaked or sperm whale) which DO create broadband pulses at regular intervals, but I've never seen them that regular.

This strongly suggests the recorded signal is man-made. The signal duration is consistently between 10 and 15 ms in duration. Highly probably it is from the Underwater Locator Beacon due to the very close match to the expected ping rate of 0.9Hz or 1 ping per 1.111 seconds, and the expected ping duration of 10ms.

However it COULD be from other repetitive signal sources such as a stray echo sounder nearby; although that would be a soul crushing disappointment to all onboard ADV Ocean Shield and the other involved in the search. Also.. WHAT CRAZY PERSON WOULD PUT A SOUNDER ON DURING AN ACOUSTIC SEARCH!?!

One last thing to note about the observed ping intervals is that they appear to increase slightly over time (order of 4ms). Whilst the sample set isn't long enough to show consistency, I see patterns like these all the time when tracking moving items fitted with pingers from stationary receivers (conceptually the same thing as the stationary pinger and moving receiver on ADV Ocean Shield).

It could be either:

- The decreasing battery voltage is causing the ping interval to increase. However it doesn't seem likely this phenomenon would be observed over such a short time window

- The audio recording segment corresponds to when the towed pinger locator has performed a Closest Point of Approach (CPA) and is now opening in range from the source. When this happens a effect similar to Doppler occurs in that the time between successive pings appears to increase as because between successive pings the receiver has moved further away.

For the segment analysed it corresponds to a slant range increase of 6 metres over 20 seconds, or 0.3m/s or 0.6 knots. If the towed pinger locator is transiting at 2 knots and the depth of the receiver is known, you could estimate a very rough range to the source by looking at the vector component in the horizontal direction if you assume the vessel passed over it. If you had the entire data set you could do much more interesting analysis!

On other forums I've had people ask about using Doppler to localise the pinger. Due to the pulse length (10ms) there is an inherent 100Hz bandwidth in the signal, which is equivalent to the frequency shift you would expect if moving at 4m/s. In other words, Doppler measurements of a single pulse is NOT going to be useful in attempting to localise.

However, because the pinger has a relatively stable repetition rate (about 1.106s) you can do a cool trick to localise it (assuming the source is stationary):
  • Perform a pass over the area of interest (ie where the pinger is)
  • Measure the ping times of arrival (either manually with waveform inspection or write an envelope detector)
  • Plot the ping times of arrival versus the time of arrival modulo ping interval. ie in MATLAB: plot(toa, mod(toa, 1.106), '.');
  • The plot indicates the relative distance from the source / receiver. Hopefully the plot should contain a local minimum representing a closest point of approach (ie like http://imgur.com/v1ZclaX). This means that the source lies in a line perpendicular to the vessels course centred at the location of the receiver (ie TPL-25) at CPA time
  • Perform another pass over the area of interest on a course perpendicular to the first
  • Measure the ping times of arrival and look for CPA as before
  • You should now have two intersecting lines. Guess what should be at the intersection? :)
Even if a consistent sequence of pings over CPA is not measured, a single segment can provide an estimate of a bearing to the source due to the reduction (towards) or increase (away) of the ping interval. This relies on having a stable ping repetition rate.

Edit: I have added a new post where I wrote a simple MATLAB script to automatically detect and visualise the MH370 pings: Automatic detection of MH370 ULB pings

Update 19th Apr - Added zoomed in plot of a single ping.

Saturday, June 19, 2010

Non-blocking Output Socket

My MATLAB TCP/IP example has been quite popular, and I myself use it at least once a week on various applications. However, one of the main limitations when compared to our C++ style output sockets was the limitation of only 1 client per socket.

What was worse, is that I knew it COULD be done. I had just not found the time to write it! Cue a spare Saturday!


The Problem

As per my previous server.m examples, you begin by constructing a Java ServerSocket.

>> server_socket = ServerSocket(output_port); 
>> server_socket.setSoTimeout(1000);

This attempts to bind the ServerSocket to the specified port. HOWEVER, it is not listening for incoming connections! To do this, you must tell the ServerSocket to accept().

>> output_socket = server_socket.accept();

This is a BLOCKING call. That is, the code will wait for either the timeout (set above) to occur, or a client to connect.

For a lot of what I do, the output socket exists mainly as a debugging / monitoring interface. So I don't even care if noone is listening. So what we need is something that manages the connection of clients and the distribution of any written data to each of them.

Enter OutputSocket and ListeningThread.

OutputSocket

OutputSocket provides a very simple interface:

//--------------------------------------------------------------------- 
// Description : Create an OutputSocket on the specified port

//---------------------------------------------------------------------

public
OutputSocket(int port) throws IOException

//---------------------------------------------------------------------
// Description : Write the supplied byte array to all connected clients

//---------------------------------------------------------------------

public void
write(byte[] data)


//---------------------------------------------------------------------

// Description : Close the OutputSocket. This closes all client

// connections and stops listening for new connections.
//---------------------------------------------------------------------

public void
close()



ListeningThread

The OutputSocket utilises ListeningThread to manage the accepting of new connections. ListeningThread implements the Runnable interface so can be started in its own thread (vital!).

//-------------------------------------------------------------------------
// Description : Helper class that listens for new TCP/IP client
// connections on a ServerSocket
//
// Parameters : server_socket - ServerSocket to listen on
// connection_list - List of connected client
// Sockets to update
// connection_stream_list - List of connected client
// DataOutputStreams to update
//---------------------------------------------------------------------
public ListeningThread(ServerSocket server_socket,
List connection_list,
List connection_stream_list)


//---------------------------------------------------------------------
// Description : Called on Thread.start()
//---------------------------------------------------------------------
public void
run()


//---------------------------------------------------------------------
// Description : Stop the thread
//---------------------------------------------------------------------
public void
stop()



Interaction

The OutputSocket constructor creates the ServerSocket as before and stores it as a member variable. It also creates 2 empty lists, one for the Socket client connections and one for the associated DataOutputStreams used to write to them.

It supplies a reference to these to the ListeningThread on construction and then starts it running in its own Thread. This is the important part. OutputSocket and ListeningThread both share a reference to the DataOutputStream list, so as ListeningThread creates connections and output streams, these are pushed onto the client list.

When OutputSocket::write() is called, it writes the data across all output streams (be there 0 or 100!)

When you have finished with the OutputSocket, a call to close() will stop the ListeningThread from running and close all client connections.

Usage from within MATLAB

Place the .class files in a suitable location (I've stuck them in my current working directory).
>> % Add the directory containing the .class files to your Java path
>> javaaddpath(pwd)
>> % Create the OutputSocket
>> output_socket = OutputSocket(1234)

output_socket =

OutputSocket@1c68b20

>> % Write some data to the output socket - text will do for now
>> for i = 1:100
output_socket.write(int8(sprintf('This is line %03d\r\n', i)));
pause(1);
end
>> % cleanup the socket - close connections
>> output_socket.close()
>> clear output_socket % no longer needed / useful


As you can see, multiple connections are supported. They can connect/disconnect at any time. Even if there are no connections, a call to OutputSocket::write() will be successful.


**Update**

I've included in the revised submission to MATLABCentral a wrapper MATLAB class that helps manage the lifetime of the Java object.

This is VERY useful because if you happen to clear the Java object (ie clear output_socket) before calling .close() the Thread and ServerSocket etc remain in memory (and stay connected!) until MATLAB is closed!

Suggested usage of MatlabOutputSocket below.
>> javaaddpath(pwd); % get the Java class
>> output_socket = MatlabOutputSocket(1234);
for i = 1:100
output_socket.write(int8(sprintf('Line %03d\r\n', i)));
pause(1);
end
>> output_socket.close()
>> clear output_socket

I've included the Java source file in the MATLABCentral submission, so feel free to have a look at the source code, have a fiddle, recompile.

To do this you will need to install the Java Development Kit (JDK).

To compile the classes, navigate to the directory containing the .java file and type:
>> javac OutputSocket.java
NOTE: You will have to repeat the call to javaaddpath() for the new files to be seen by MATLAB.

Monday, September 7, 2009

TCP/IP Socket Communications in MATLAB - Part 2

So my little blog post about TCP/IP Socket Communications in MATLAB was rather popular despite being a very simple example.

One of the main limitations that people found when trying to utilise the server/client scripts for their own applications, was that it was incredibly inefficient at shifting large volumes of data around. This is thanks to the following lines:
message = zeros(1, bytes_available, 'uint8');
for i = 1:bytes_available
message(i) = d_input_stream.
readByte;
end
For each individual byte within the message, there was the overhead of a function call.

I mentioned in comments on MATLAB Central that I made a Java class that bypassed this overhead and allowed for more efficient transfer. Being in a caring sharing mood (and getting sick of emailing it to people all the time!) I thought I would upload it also.

The is available on the Mathworks File Exchange: TCP/IP Socket Comms Example using Java Class

Compiled Java?

The client/server scripts are essentially identical to their previous iterations.

The server still uses a ServerSocket and the resulting DataOutputStream to which we write data.

The client side uses a Socket to connect to the specified host and port which provides us an InputStream which we wrap in a DataInputStream to read data from. However, instead of using the interface of DataInputStream directly, we hand off to another function (the new DataReader class) to perform the read.

The code for the example client is outlined below. I have bolded the changes. (The server is unchanged).

client.m
% CLIENT connect to a server and read a message
%
% Usage - message = client(host, port, number_of_retries)
function message = client(host, port, number_of_retries)

import java.net.Socket
import java.io.*

if (
nargin < number_of_retries =" 20;" style="color: rgb(0, 153, 0);">% set to -1 for infinite
end

retry = 0;
input_socket = [];
message = [];

while true

retry = retry + 1;
if ((number_of_retries > 0) && (retry > number_of_retries))
fprintf(1, 'Too many retries\n');
break;
end

try
fprintf(1, 'Retry %d connecting to %s:%d\n', ...
retry, host, port);

% throws if unable to connect
input_socket = Socket(host, port);

% get a buffered data input stream from the socket
input_stream = input_socket.
getInputStream;
d_input_stream =
DataInputStream(input_stream);

fprintf(1, 'Connected to server\n');

% read data from the socket - wait a short time first
pause(0.5);
bytes_available = input_stream.available;
fprintf(1, 'Reading %d bytes\n', bytes_available);

data_reader = DataReader(d_input_stream);
message = data_reader.readBuffer(bytes_available);


message = char(message'); % Data comes out as a column vector

% cleanup
input_socket.close;
break;

catch
if ~
isempty(input_socket)
input_socket.close;
end

% pause before retrying
pause(1);
end
end
end

Instead of looping for each byte, we now ask the DataReader object to read the specified number of bytes and return them to us. So a single function call per read. This could be used to read buffers of expected data sizes.

DataReader

Below is the Java source for the DataReader class (available in the MATLAB File Exchange also).

import java.io.*;

class DataReader
{
public DataReader(DataInput data_input)
{
m_data_input = data_input;
}

public byte[] readBuffer(int length)
{
byte[] buffer = new byte[length];

try
{
m_data_input.readFully(buffer, 0, length);
}

catch (StreamCorruptedException e)
{
System.out.println("Stream Corrupted Exception Occured");
buffer = new byte[0];
}
catch (EOFException e)
{
System.out.println("EOF Reached");
buffer = new byte[0];
}
catch (IOException e)
{
System.out.println("IO Exception Occured");
buffer = new byte[0];
}

return buffer;
}

private DataInput m_data_input;
}
It is very simple class that is given a DataInput to read from, and has a single public method readBuffer(int length) that returns a byte array.

If you have the Java Software Development Kit installed, you can compile the class with:
C:\matlab\matlab_socket> javac data_reader.java

Running the Client/Server example with DataReader

This is the same as the previous example, with one exception. You must tell MATLAB where to find the compiled DataReader class (In my case it is located in the C:\matlab\matlab_socket directory):
>> javaaddpath('C:\matlab\matlab_socket');
Opening up two instances of Matlab:
% Instance 1
>> message = char(mod(1:1000, 255)+1);
>> server(message, 3000, 10)
Try 1 waiting for client to connect to this host on port : 3000
Try 2 waiting for client to connect to this host on port : 3000
Try 3 waiting for client to connect to this host on port : 3000
Try 4 waiting for client to connect to this host on port : 3000
Client connected
Writing 1000 bytes

% Instance 2 (simultaneously)
% NOTE: If the 'server' was runnning on a non local machine, substitute its IP address
% or host name here:
% data = client('10.61.1.200', 2666); % To connect to server at IP 10.61.1.200:2666
>> javaaddpath('C:\matlab\matlab_socket');
>> data = client('
localhost', 3000)
Retry 1 connecting to
localhost:3000
Connected to server
Reading 1000 bytes

data =



 !"#$%&'()*+,-./0123456789:;<=>?@ABCDEFGHIJKLMNOPQRSTUVWXYZ[\]^_`abcdefghijklmnopqrstuvwxyz{|}~     

 !"#$%&'()*+,-./0123456789:;<=>?@ABCDEFGHIJKLMNOPQRSTUVWXYZ[\]^_`abcdefghijklmnopqrstuvwxyz{|}~     

 !"#$%&'()*+,-./0123456789:;<=>?@ABCDEFGHIJKLMNOPQRSTUVWXYZ[\]^_`abcdefghijklmnopqrstuvwxyz{|}~     

 !"#$%&'()*+,-./0123456789:;<=>?@ABCDEFGHIJKLMNOPQRSTUVWXYZ[\]^_`abcdefghijklmnopqrstuvwxyz{|}~    

Why not call readFully() from within client.m?

The more astute readers may be asking the above question. Why not do something like this in client.m:
message = zeros(1, bytes_available, 'uint8');
d_input_stream.
readFully(message, 0 bytes_available);
At first glance that seems great. It runs without error, but hey, you don't seem to be getting any output?¿?

The problem is that DataInputStream.readFully() takes a reference to a byte array to populate.

When you pass things around in MATLAB that are to be modified, MATLAB supplies a copy. Thus the Java most likely sees a byte array reference, and populates it, but it is only the copy that is populated, not the original message array.

Until there exists some way in MATLAB to pass references to Java methods, I am stuck using helper Java classes (If anyone else has ideas on how this could be done better, feel free to comment!).

Thursday, September 18, 2008

Mathworks File Exchange Profile Information - Part 2

Ok, A quick little update here before I head off on holidays!

Continuing on from my previous post I realised that I wanted to watch how my submissions were performing over time (my memory is not quite what it used to be).

I had a crazy idea, could I create a Timer that would scrape my user statistics from the Mathworks File Exchange site and save them in a format such that I could plot them easily???

Introducing Monitor File Exchange Statistics.

Plotting my stats for the past week:
I'm not a big fan of the look of the submitted files plot as it can get a bit busy. I'm envisaging a lightweight GUI with a drop-down selector for each submitted file might be nicer. Particularly for those people with a large number of submitted files. I'll update the submission after I get back to do something like this.

To set the automatic scraping of your File Exchange statistics you need to do the following:
  1. Download the required MATLAB files from here. Need Monitor File Exchange Stats and Retrieve File Exchange Profile Information.
  2. Place the files somewhere in your MATLAB path. Note, if you already have a startup.m, just add the InitialiseFileExchangeInfoTimer and UpdateFileExchangeInfoTimerFcn to your existing startup.m and call InitialiseFileExchangeInfoTimer from within.
  3. Modify FEX_AUTHOR_ID and UPDATE_PERIOD with your File Exchange Author ID and your desired scrape frequency (seconds). The cache file defaults to being stored in your MATLAB root directory, modify FEX_SAVE_FILE if you wish to customise this location.
  4. Save any changes and restart MATLAB.
  5. Each time you restart Matlab you will be presented with a plot of your current and historical stats (for as long as you have been scraping).
  6. To manually plot your current stats, call plotLocalFileExchangeInfo(cache_file)
In terms of file size, I've had the script running every 10 minutes for nearly a week and my cache file is only 45 kB.