Spectrograms and Whale Sounds

Part 1: Build Your Own Spectrogram

Whales are very expressive, and their calls are easiest to visualize over time and frequency together. To see a signal plotted in time and frequency together, let's build a spectrogram, which is sometimes also called a short-time Fourier transform (STFT). The spectrogram copies chunks of the time signal, and does a Fourier transform on each chunk. In that way, you can see how the frequency content of the signal is changing over time. Here's what the spectrogram function needs to do:

• (a) Define a time window of length w. This magnitude of this window is 1 for the first w samples and 0 elsewhere: W(n) = 1 for n in the range [0,w], and W(n)=0 otherwise.
• (b) Multiply W(n) by input signal x, which is length L. In general L is much larger than w.
• (c) Take the Fourier transform of the result of step (b). We really only need to perform a transform of length w, not length L, since only w samples are non-zero.
• (d) Store the magnitude from step (c) in the first column of a matrix X. Let's assume that our input signal is real. This means that the magnitude of the Fourier transform is symmetric and thus, we only have to store (w/2)+1 components to fully specify the Fourier Transform of x.
• (e) Move the window over by m samples. Then we can define our new window W(n,m) = 1 for n in [m+1,m+w] and 0 otherwise. As you can see, the original W(n) was a special case of W(n,m) with m = 0. We can define o = w - m as the amount of overlap in successive FTs. This will be important to consider in our practical implementation of a spectrogram function.
• (f) Repeat steps (b) through (e), each time storing the result in the next column of X.

Create a new M-file in MATLAB. The first line of a function should define the input vectors, the output vectors and the function name:

function X = spectrogram(x,w,o,fs)

X will be the output matrix, x is our input vector, w is the length of the window, o is the length of the time overlap and fs is the sampling frequency. (Note that MATLAB has its own spectrogram command, which you will be replacing with your function. The MATLAB function has more options than you will explore here.)

The spectrogram only needs a few lines of code:

counter = 0; for k = 1:w-o:length(x)-w     counter = counter + 1;     temp = abs(fft(x(k:k+w-1)));     X(:,counter) = temp(1:w/2+1); end

Comment each line of code to explain what it does. Explain in your comment for the last line of the for-loop why we keep only half of the magnitude of the Fourier transform.

Here is plotting code. Comment each line of this code to explain what it does.

T = 1/fs; yvec = 0:(fs/2)/(w/2):(fs/2); xvec = 0:T*(w-o):T*(length(x)-1);  figure() imagesc(xvec,yvec,X) axis xy xlabel('time (s)') ylabel('frequency (hz)')

Part 2: Whale Calls

Load the calls of the two whales from whale_calls.mat (type load 'whale_calls.mat'). This file contains three MATLAB variables: X1, X2 and fs. The variable fs is the integer sampling frequency in Hertz. The variables X1 and X2 are matrices where each row is a whale call. Listen to a few from each matrix. Describe any differences you hear between the whales from X1 and X2.

Let's add some code to plot the whale songs.

w = 256; o = 128; x1 = X1(1,:); x2 = X2(1,:); spec1 = spectrogram(x1,w,o,fs); spec2 = spectrogram(x2,w,o,fs);

Comment each line of code above.

View different whale calls from X1 and X2. Be ready to explain how what you see corresponds to what you hear. Play around with w and o and describe what changing each one does to the quality of the resulting image. What happens when you set w = o? Why? Be ready to explain the effect of variables w and o.

Part 3: Multipath in the Ocean

Echoes occur not just in canyons, but can be a challenging part of any remote-sensing signal processing system, including radar, ultrasound, and sonar. When the transmitter and receiver are at different locations, the analog of echoes is that the receiver hears the transmitter's signal, and then hears echoes caused by the signal traveling different (usually longer) paths. This is called multipath, and it can be well-modeled by a linear time-invariant filter.

Imagine that you are on a submarine with sonar sensors. Somewhere in the distance, a whale emits a call and your microphone picks it up. What would your recording of the whale sound like? The call of the whale goes through multipath distortion. Let us assume that we have a two-path multipath. Since the whale's vocal cords create an omnidirectional signal, the microphone picks up the signal that is on a linear path from the whale, (called line of sight, for obvious reasons), and picks up a second version of the signal that first bounces off the ocean floor. If we call the line of sight signal x(t), then the signal that bounced, received T seconds later, would be a*x(t-T) (where "*" means times). The constant 'a' (between 1 and -1) accounts for the attenuation that occurs when the signal loses energy by hitting the ocean floor (it can be negative due to the possibility that the reflection turns the signal upside down). The delay term 'T' is needed because the path taken is necessarily longer than the linear one. In practice, there may be many bounces of the signal off the sea floor and off the ocean-air interface (and possibly off other whales or sunken ships, etc.). So, a general expression for the received signal would be:

Figure 1

An even more general expression would add a weight a0 in front of the x(t) to allow some attenuation of the line-of-site term.

Does that formula look familiar? It should... It's convolution! We can represent the delays and attenuations with an impulse response:

Figure 2

and the received signal z(t) as h(t) convolved with x(t). In MATLAB, we will represent this kind of continuous-time impulse response by sampling it with sampling interval Ts. Then the delay times T1..Tn are represented as multiples of the sampling interval, so for M samples:

Figure 3

Many of the ak's will be 0 if there are no delay components corresponding to (k)Ts seconds. To specify the impulse response, we need only specify the ak values for k=0...M (assuming we know Ts). Here's an example of a multipath h(t) : [1 0 0 0 0 .6 0 0 -.45 0 0 .1]. The numbers in the list represent the weights [a0 a1 . . . aM] where in this case M=12.

Figure 4: A ship is shown with a sensor array descending from it, an object is denoted in red (whale? dolphin? submarine?) and the object emits signal x(t), which undergoes multipath. The total effect of the multipath is modeled as a channel with impulse response h(t). The ship sonar system receives z(t) = x (t) * h(t). (In practice, the ship really receives z(t) = x(t) * h(t) +n(t), where n(t) is some extra noise.)

Multipath in the Ocean

Figure 5: This is an example of the impulse response for the multipath in a 50 meter water channel (for example, in a bay).

Impulse Response Example

Can you think of a case when a0 would be 0? Hint: in this case there would be no direct path between the receiver and the transmitter. Can you explain why most of the weights ai in the list for h(t) are zero (why h(t) is sparse)?

Now we're going to simulate multipath and see how it changes the whale calls. To get a sparse h(t), we'll use a random number generator (rand, which generates numbers uniformly over [0,1]) and scale the shifted [-0.5,0.5] version so that only rarely will the round function give 1 or -1. In the code below, the closer that sparseness is to 1, the more sparse the impulse response will be.

sparseness = 1.0005; scaling = .2-.00001:-.00001:0; h = [1 scaling.*round(sparseness*(rand(1,20000)-.5))];

Plot the impulse response of this filter and comment the lines of code.

Create a new impulse filter that has exponentially decreasing coefficients as opposed to linear ones. (Hint: you only have to modify the scaling vector). Plot it.

Create z1 = conv(h,x1) and z2 = conv(h,x2). Listen to these two signals and visualize them with both time and spectrogram plots. How did the multipath affect the whale calls? Play around with the formula and describe how changing various parts of it affects the distortion. Can you find a range for the sparseness coefficient such that you can no longer tell which whale is which?