# Graphical representation of a moving sound source

Some explanation since I have got a specific and another more open-ended question: I used two microphones in a 90° arrangement (see picture below) to capture sound from a source moving around the room in a circular /elliptical fashion. The intention is to get an idea about directional hearing. (A second measurement with a KEMAR dummy-head was done as well).

I saved the measured data as .wav (48kHz) and imported the .wav file into Mathematica. (The .wav- file can be found here). That made a pretty long list - approximately 900000 entries per channel. Hence I used

 dropfunction[data_, n_] :=
Table[Nest[Drop[#, {1, -1, 2}] &, data[[i]], n], {i, 1, 2}]


to reduce the resolution. I found by hearing/looking at ListLinePlot that n = 7 is a good catch. That leaves around 7000 entries per channel.

Here comes my first and very specific question: I plotted the absolute values of my list in a fashion like:

ListLinePlot[Abs@data[[1]], -Abs[data[2]]]


and got:

with left and right channel clearly divided along the x-axis.

Question 1: Is there a possibility to flip the plot in such a way that x-axis becomes the y-axis? (So that for instance the "loudness" of the right channel is on the right hand side). I tried PairedBarChart actually gives the kind of representation I wish, but It is just not suitable for my data.

Question 2: Does anybody else have a great idea how to graphically represent my data? I tried plotting the difference of the absolute values of both channels to get a graph showing which channel is currently louder at the positions in time/on my elliptical path. With some lowpass-filtering I got:

-
Just a note: Part may help you thinning your list. E.g. x[[1;;-1;;10]] takes every tenth element from head to tail. –  Yves Klett Mar 19 '13 at 17:28
@Sascha As a general note, ListLinePlot returns a regular Graphics object. There isn't anything special about the "built-in" plot functions (as far as the Mathematica system is concerned) and you can build your own special plots pretty easily. –  amr Mar 19 '13 at 17:29
As for plotting, perhaps a ListPolarPlot will be enlightening. –  Yves Klett Mar 19 '13 at 17:31
I created my downsampled data by using the posted function dropfunction (with n= 7) on the uploaded .wav-file (or as Yves Klett suggested by using Part) –  Sascha Mar 19 '13 at 19:37
V9 has the Downsample function. –  Sjoerd C. de Vries Mar 20 '13 at 6:22

To visualize directional information between the two signals, it is sufficient to cross-correlate them and look at the time lags. If the signal in, say, microphone A lags that in microphone B, then it implies that the source was closer to B than A. As you move around the room, the lag should change appropriately depending on the sine/cosine (depending on the convention) of the angle subtended between line of the microphones and the source.

So moving on to processing, I'll first low pass filter and downsample your data, since most of the frequency content is in the lower frequencies. Note that simply downsampling (as in your nested Drop) without filtering will result in higher frequencies aliasing into the lower ones.

{fs, {ph1, ph2}} = Import["~/Downloads/T55test.WAV", #] & /@ {"SampleRate", "Data"};

Clear@decimate
decimate[data_, fs_, n_] := {fs/n, Downsample[LowpassFilter[data, fs/n π, SampleRate -> fs], n]}

{dfs, dph1} = decimate[ph1, fs, 6];
{dfs, dph2} = decimate[ph2, fs, 6];


Here I've downsampled by a factor of 6, meaning the new sampling frequency is 8000 Hz. Next, I chop the data into chunks of 5 milliseconds with 50% overlap between consecutive chunks and cross-correlate them.

Clear@corr
corr[{d1_, d2_}] := ListCorrelate[d1, d2, {-1, 1}, 0]

With[{seg = dfs/200}, With[{h = N@Array[HammingWindow, seg, {-1/2, 1/2}]},
c = Composition[Rescale, corr] /@
Transpose[DeveloperPartitionMap[h # &, #, seg, seg/2] & /@ {dph1, dph2}];
MatrixPlot[Transpose@c, ColorFunction -> "Temperature",
FrameTicks -> {{-1000. seg/dfs/2, 0, 1000. seg/dfs/2}, {0, 5, 10, 15}, False, False},
AspectRatio -> 1/5, FrameLabel -> {"Lag (ms)", "Time (s)"}, DataReversed -> True,
DataRange -> {{0, Length@dph1/dfs // N}, {-1000. seg/dfs/2, 1000. seg/dfs/2}}
]
]]


You can now clearly see the lag — it is zero at the start and end, when the source is symmetrically positioned between the two phones and increases and decreases as the source moves. You can now make out the direction from this plot, or perhaps implement a simple delay-and-sum beamforming to get the direction. Note that you can also do the same processing with the original data, change overlap (or none at all), change the windowing, etc., and the basic structure of the lag will remain the same.

As an interesting exercise, one can also make a guess as to what the spacing between the two microphones might be from looking this lag plot. Since the speed of sound in air is approximately 330 m/s and by eyeballing the plot and approximating the max lag to be at around 1 ms, one can infer that the two microphones must be separated by roughly 30 centimeters (1 ft).

-
I think this is an excellent answer. What is missing is the definition of the correlation matrix. –  Stefan Mar 20 '13 at 10:01
What Stefan means is that one fails to fine a definition of corr when looking at your code. –  Sascha Mar 20 '13 at 21:00
@Sascha Oh crap! Didn't notice that I missed that. I've added it now :) –  The Toad Mar 20 '13 at 21:11
@Sascha btw, is that estimate for the separation between the microphones correct? –  The Toad Mar 20 '13 at 22:03
Ah yes, those are highly directional, so interpreting the lag is a lot more complicated and will also need to take into consideration the room geometry... –  The Toad Mar 22 '13 at 14:56

How about plotting it around a circle? This just wraps the amplitude signal on top of a circle which represents the time of your signal.

data = Import["directory/T55test.WAV"][[1, 1]];
signal = Transpose[Norm /@ Partition[#, 480] & /@ data];
s1 = Transpose[signal][[1]];
s2 = Transpose[signal][[2]];
aroundcircle[datain_] := {(1 + 10 datain[[#]]) Sin[(# 2 \[Pi])/
1894], (1 + 10 datain[[#]]) Cos[(# 2 \[Pi])/1894]} & /@ Range[1894];
s1around = aroundcircle[s1];
s2around = aroundcircle[s2];
Show[ListLinePlot[{s1around, s2around}, PlotStyle -> {Red, Blue}],
Graphics[Circle[{0, 0}, 1]], AspectRatio -> 1,PlotRange -> {{-1.5, 1.5}, {-1.5, 1.5}}]


The value 1894 is the number of elements in the sampled signal.

-
NB. I've taken the second line from @Simon Woods' code. You can take any filter or sampling of the data and the plot method will be the same. –  Jonathan Shock Mar 20 '13 at 3:18
You can also define the ellipse precisely as it is rather than the simple circle. –  Jonathan Shock Mar 20 '13 at 3:20
That is a really nice idea:) –  Sascha Mar 20 '13 at 8:58
I should also note that in this code, nothing about the 90 degree offset of the microphones is encoded. This reduces the usefullness of the code but it could still be adapted. –  Jonathan Shock Mar 20 '13 at 9:25
If I run your code on my machine, your result isn't reproduced, any idea why that could be? –  Sascha Mar 20 '13 at 9:31

Let's begin by considering the purpose of the visualization. It is primarily to assess direction of the sound source. The most useful visualizations therefore will

1. Provide ways to assess relative amplitudes (rather than differences)--they will likely use logarithmic scales for that--and

2. They will clearly distinguish when the left hand microphone's response is larger than the right hand microphone's response. Because the source apparently ended where it started, circular or periodic representations of these relationships may be useful.

Now we apply visualization principles to design an effective graphic.

To enable rapid, accurate comparison of two signals $X$ and $Y$ (which, by virtue of #1, should be logarithms of the amplitudes in the two channels), we should use graphical elements that reflect the differences $X-Y$. We wish to see how they (a) relate to the collective "strength" $(X+Y)/2$ and (b) evolve over to time (as a secondary indicator). A scatterplot of $X-Y$ (log of the signal ratio) against $(X+Y)/2$ (log of the geometric mean) best achieves (a), by (i) representing both coordinates as distances along orthogonal axes and (ii) using a horizontal line as a reference; namely, the line $X-Y=0$: vertical deviations from that line will let us assess the relative strengths of the two microphones. This leaves a host of secondary symbols, such as line thickness, dashing, and color, to represent the temporal evolution. Due to the circular or quasi-periodic nature of this experiment, a periodic color scale seems suitable.

Some subtle alternation in the plot to distinguish positive values of $X-Y$ from negative values will help further. Below, I chose simply to vary the lightness of the background, wishing not to complicate the graphic any further: it already shows a lot (total signal strength, relative strength, and time for thousands of observations).

Some statistical insights offer potential improvements.

Data can be noisy. To help us see the big patterns, some amount of smoothing is advisable. A convolution with a Gaussian kernel is a quick and easy way to perform a smooth when the signals are not otherwise impaired by extreme outliers. A quick check of the data indicates there are no outliers evident.

Here are examples motivated by these design considerations, showing the results for a range of smoothing strengths.

data = Import["F:/temp/T55test.WAV"][[1, 1]];
signal = Transpose[Norm /@ Partition[#, 480] & /@ data];

Table[
kernel = Table[Exp[-i^2], {i, -3, 3, 6/n}]; (* Gaussian kernel of length around 2n+1 *)
kernel = kernel / Total[kernel];            (* Normalization to unity *)
smooth = ListConvolve[kernel, #~Join~#[[;; Floor[n]]]] & /@ (signal\[Transpose]);
ListLogLogPlot[{Sqrt[First[smooth] Last[smooth]], First[smooth]/Last[smooth]}\[Transpose],
Joined -> True,
ColorFunction -> (Hue[2/3 #2] &),
PlotStyle -> Thickness[0.01],
PlotLabel -> "Smooth: " <> ToString[n],
AxesLabel -> {"Geometric\nmean", "Ratio"},
Prolog -> {GrayLevel[.95], Polygon[Log[{{0.001, 1}, {0.1, 1}, {0.1, 0.1}, {0.001, 0.1}}]]}],
{n, {1, 301, 1001}}]


(The visual metaphors at play in these graphs work nicely, because there is a strong association of time with relative signal strength: cool colors (blues and greens) accent the points where $X$ exceeds $Y$ and warm colors (reds, oranges, and yellows) indicate where $Y$ exceeds $X$. Thus, color, height above or below the visual horizon, and darkness of background all cooperate in helping us see the aspects of these data that are most important to the analytical objective.)

The raw signal on the left appears to bear little relationship to the path of the sound source. However, the heavily smoothed signal on the right is promising, showing a strong qualitative relationship to the location of the sound. The middle image preserves some of the major details.

To go further, consider replacing Table by Manipulate to explore the effects of the smoothing.

-

Following on from @Simon Woods, adding a colour function to the plot allows you to see the data evolve with time (e.g. from yellow to red):

    data = Import["T55test.WAV"][[1, 1]];
signal = Transpose[Norm /@ Partition[#, 480] & /@ data];

tpos = Range [Length[signal]];

tcolour = Blend[{{Min[tpos], Yellow}, {Max[tpos], Red}}, #] &;
tcolplot = Graphics[MapThread[{tcolour[#1], Point[#2]} &, {tpos, signal}],
Axes -> True, AspectRatio -> 1]


-

How about plotting the signal strength from one microphone on the $x$ axis, and from the other on the $y$ axis.

data = Import["T55test.WAV"][[1, 1]];
signal = Transpose[Norm /@ Partition[#, 480] & /@ data];
ListPlot[signal, AspectRatio -> Automatic]


-

Maybe just use ListPlot instead of ListLinePlot as follows:

dataLeft = Table[Sin[1000 x] + Sin[1005 x], {x, 0, 1000}];

dataRight = RotateRight[dataLeft, 100];

ListPlot[{Transpose[{#, Range[Length[#]]} &@(-Abs /@ dataLeft)],
Transpose[{#, Range[Length[#]]} &@(Abs /@ dataRight)]},
Joined -> True]


-
Do you have any idea how to change the Ticks or FrameTicksso that they show only a absolute distance from 0 instead of +/- ? –  Sascha Mar 19 '13 at 18:22
Easiest way could be specifying custom ticks: ListPlot[{Transpose[{#, Range[Length[#]]} &@(-Abs /@ dataLeft)], Transpose[{#, Range[Length[#]]} &@(Abs /@ dataRight)]}, Joined -> True, Ticks -> {{{-2, 2}, {-1, 1}, 1, 2}, Automatic}] –  halmir Mar 19 '13 at 19:25
@halmir Thanks - that's what I'd do too. –  Jens Mar 19 '13 at 19:53

Combining Simon's and Jonathan's ideas with a small twist:

data = Import["T55test.WAV"][[1, 1]];
signal = Transpose[Norm /@ Partition[#, 480] & /@ data];

{intf1,  intf2} = Interpolation[10 signal[[All, #]]] & /@ {1, 2};

pplt1 = ParametricPlot3D[{
{ 3/2 Cos[(t 2 \[Pi])/1894], Sin[(t 2 \[Pi])/1894],  intf1[t]},
{ 3/2 Cos[(t 2 \[Pi])/1894], Sin[(t 2 \[Pi])/1894], intf2[t]}},
{t, 1, 1894}, BoxRatios -> 1,  Mesh -> None, ColorFunction -> (Hue[#3] &),
PlotStyle -> {Thick, Directive[Thick, Dashed]},
PlotRange -> {{-2, 2}, {-1, 1}, {0, .6}}];

pplt2 = ParametricPlot3D[{3/2 Cos[(t 2 \[Pi])/1894], Sin[(t 2 \[Pi])/1894],  v},
{t, 1, 1894}, {v, 0, 1},
BoxRatios -> 1, PlotStyle -> Opacity[.1], ColorFunction -> (Hue[#3] &), Mesh -> None,
PlotRange -> {{-2, 2}, {-1, 1}, {0, .6}}, ImageSize -> 600];

Show[pplt2, pplt1]
`

-
Looking at your graphic I came to think it would be awesome for actual presentation to print a banner with the data in 2D and then wrap it up to get a cylinder barrel of paper much like your plot. You could then stick your head in the ring and view the data in relation to its actual physical position. Sad that I don't have to present it though (it's only a report to be handed in) :) –  Sascha Mar 21 '13 at 11:13
Sacha, if you are allowed to turn in your report in CDF form, you actually can get pretty close to "sticking your --or, better yet, your teacher's -- head in the cylinder barrel":) In case you haven't already seen it, you should watch Vitaly Kaurov's presentation Mastering Dynamic Presentations to see how that can be done. –  kguler Mar 21 '13 at 11:28
unfortunately they want a .pdf-file and I suspect just for the sake of printing it on top of that. (And I currently don't know if the hassle of getting 3D data from Mathematica to .pdf is worth the hassle anyway, meaning that I suspect there to be some kind of quality issue, color, aliasing etc.) –  Sascha Mar 21 '13 at 11:34