Sound speed measurement (with two microphones)

I have two microphones connected to my laptop through one microphone port. The microphones are separated by some distance (about 1.5 meters). My aim is to measure the delay of sound between microphones to measure sound speed. (v_sound=x/delay). Below is my program:

if true
 clear; figure; grid on; hold on; 
Fs = 4000; % sampling rate in Hz 
dt = 8;    % duration in seconds 
% get data 
y = wavrecord(dt*Fs, Fs, 2); 
% convert to time (sec) 
tmax = length(y)/Fs; 
t = linspace(0, tmax, dt*Fs); 
% plot 
plot(t*1000,y); 
axis([0 tmax*1000 -2 2]); 
xlabel('time (ms)'); 
% --- Analysis --- 
j = 1; 
for i=1:length(y) 
if y(i)>0.15 
    pick(j) = 1000*i/Fs; 
    fprintf('%3d --> %8.1f ms\n',j, pick(j)); 
    j=j+1; 
end 
end
end

This program runs 8 seconds, during that time i did two claps with my hands close to one microphone (so one microphone receives sound before the other, i need that delay). And here is an output:

Blue color is a microphone close to my hands and green is a microphone that is abou 1.5m. away. I zoom in and try to find the delay manually using data cursor but there are too many peaks there, i don't know which pair of peaks to take.

Do you have any suggestions? Ideas?

[Star Strider]

Post your data (use the ‘paperclip’ icon to choose and attach the file. We can then experiment with it and suggest solutions.

Since the time is relative, uploading the segment between 2000 ms and 5000 ms would be best. By my calculations, the difference between the microphone recordings of the same impulse would be only about 4 ms or 17 samples.

[Akbar]

yes your calculations are correct. Before uploading this file i did a two clap experiment again, just so you know. And the exact distance between the mics is 1.66 meters.

[Star Strider]

Thank you! I can get the data from your .fig file.

[Akbar]

it would be great if MATLAB will calculate delay for me. Below Image Analyst wrote some piece of code. Could you help me with its implementation to code that is above and i will try running it.

Answer 1

[Star Strider]

Success!

It took a few minutes to get the data out of the figure, and then some experimenting. I thought about low-pass filtering your signal (I still believe that would be a good idea if you intend to repeat this experiment), but decided in the interests of time to go with simple thresholding (while I watched Manchester City and West Ham).

I include all the code, including that to acquire the data from your .fig file (even though you don’t need that part).

The code:

% Get figure and data, convert to numeric from cell:
ssf = openfig('sound_speed.fig')

hf2 = get(gca, 'Children')
hf3xc = get(hf2, 'XData')
hf3yc = get(hf2, 'YData')

hf3xd = cell2mat(hf3xc);
hf3yd = cell2mat(hf3yc);

trshld  = -0.6;                  % Threshold arbitrary, estimated from plot

% Threshold first:
pulse1 = find(hf3yd(1,:) < trshld);
pulse2 = find(hf3yd(2,:) < trshld);

Fs = 4000;                       % Sampling frequency (Hz)

% Equalise vector lengths: 
pl1 = length(pulse1);
pl2 = length(pulse2);
vlen = min(pl1, pl2);
pulsx = pulse1(1:vlen);         % Lengths in samples
pulsy = pulse2(1:vlen);
pulsxt = pulsx/Fs;              % Lengths in seconds
pulsyt = pulsy/Fs;

% Subtract:
PulsDs = pulsx - pulsy;         % Difference in microphones arrival times (samples)
PulsDt = pulsxt - pulsyt;       % Difference in microphones arrival times (seconds)

% Take mean:
PulsDsM = mean(abs(PulsDs));
PulsDtM = mean(abs(PulsDt));

MicDist = 1.66;                 % Distance between microphones (m)
Csc = MicDist/PulsDtM           % Calculated speed of sound (m/s)

fprintf(1,'\n\n\tSpeed of sound  = %.3f m/s\n\n', Csc)

The result:

      Speed of sound  = 390.588 m/s

If you intend to repeat this experiment over time, I suggest:

• Develop an appropriate low-pass filter
• Develop an adaptive threshold or other peak- or pulse-detection algorithm
• Record everything that can affect the air density and therefore the speed of sound (temperature °K, barometric pressure, relative humidity, etc.)
• Collect data over time and do regressions (and plots if the function space permits) of the various variables and how they affect the speed of sound over time
• This sounds like a cool advanced secondary or undergraduate university physics experiment!

[Akbar]

WOW! this seems to be very good answer. But i couldn't get it work. I run your program with another figure but it didn't work. Output: Speed of sound = NaN m/s. Here is a figure i've tried:

and why is your threshold -0.6 ? How did you choose this value?

[Akbar]

For another experiment i made your program to read Figure2: output is: Speed of sound = 1.724 m/s

[Akbar]

Speed of sound = 474.286 m/s

[Star Strider]

Success!

Filtering and adaptive thresholding worked. Here is the revised section of my earlier code (note the leading and trailing comments so you know where it goes with respect to the previous code). What isn’t here, didn’t change:

% Get figure and data, convert to numeric from cell:
figfile = 'Figure3.fig';
ssf = openfig(figfile);

hf2 = get(gca, 'Children');
hf3xc = get(hf2, 'XData');
hf3yc = get(hf2, 'YData');

hf3xd = cell2mat(hf3xc);
hf3yd = cell2mat(hf3yc);

% Bandpass filter:
[b,a] = butter(3, [0.2  .9])
y1f = filtfilt(b, a, hf3yd(1,:));
y2f = filtfilt(b, a, hf3yd(2,:));

pulsidx1 = find(y1f >= max(y1f)*0.5);
pulsidx2 = find(y2f >= max(y2f)*0.5);

pulse1 = pulsidx1;
pulse2 = pulsidx2;

Fs = 4000;                       % Sampling frequency (Hz)
Ts = 1/Fs;                       % Sampling interval (seconds)

% Equalise vector lengths:

The fprintf statement changed to:

fprintf(1,'\n\n\tSpeed of sound for %s = %.3f m/s\n\n', figfile, Csc)

The Butterworth filter coefficients are:

b =

     374.4527e-003     0.0000e+000    -1.1234e+000     0.0000e+000     1.1234e+000     0.0000e+000  -374.4527e-003

a =

       1.0000e+000   633.3601e-003  -999.0030e-003  -304.9157e-003   674.1771e-003    97.3976e-003  -137.7613e-003

The results I got:

  Speed of sound for sound_speed.fig = 442.667 m/s
  Speed of sound for Figure1.fig = 349.474 m/s
  Speed of sound for Figure2.fig = 339.270 m/s
  Speed of sound for Figure3.fig = 349.474 m/s

I cannot guarantee that it’s completely robust to all possible data, but it works for what we presented it with thus far.

[Akbar]

yeah, unfortunately it's not robust to all possible data.

Speed of sound for Figure5.fig = 2.473 m/s

Speed of sound for Figure6.fig = 450.169 m/s

Speed of sound for Figure7.fig = 173.974 m/s

Speed of sound for Figure8.fig = 2.150 m/s

Can we make it work for these data as well, thank you very very much, you did a lot for me!

[Star Strider]

Change the thresholds. This works for everything you’ve sent me thus far:

pulsidx1 = find(y1f >= max(y1f)*0.95);
pulsidx2 = find(y2f >= max(y2f)*0.95);

Remember, the algorithm assumes that you are directly on the axis of the microphones, not at an angle to them. I didn’t see echoes, but there does seem to be significant background noise, even after filtering.

  Speed of sound for sound_speed.fig = 415.000 m/s
  Speed of sound for Figure1.fig = 349.474 m/s
  Speed of sound for Figure2.fig = 442.667 m/s
  Speed of sound for Figure3.fig = 316.190 m/s
    Speed of sound for Figure5.fig = 368.889 m/s

My pleasure!

The sincerest expression of appreciation on MATLAB Answers is to accept the answer that most closely solves your problem. It also helps others with similar questions find applicable answers.

[Star Strider]

Thank you!

[Hamza Ashraf]

i want to do the same thing but my problem is i have recorded sound as wav format how can i find speed of that sound at differnt mircophones seperated by a distance. please help me out

Answer 2

[Image Analyst]

Can't you just threshold?

sound1TimeIndex = find(abs(sound1) > 0.2, 1, 'first');
sound2TimeIndex = find(abs(sound2) > 0.2, 1, 'first');
deltaTimeIndex = abs(sound1TimeIndex  - sound2TimeIndex);

Did you try that?

[Akbar]

no i haven't tried, how can i implement you piece of code here, let's try together.

[Akbar]

I think using this code my program will calculate delay for me, that would be great!!! But i can't implement it here, could you please help me with implementation and i will try

[Image Analyst]

Well how did you plot it? sound1 is whatever you plotted in green and sound2 is whatever signal you plotted in blue. You do have two signals don't you?

[Akbar]

yes, there are two signals. both of signals are in this line of code from my code above: y = wavrecord(dt*Fs, Fs, 2); And: plot(t*1000,y); % gives plot of both signals

[Image Analyst]

t is your time axis, which you multiplied by 1000 for some reason. y is the sound signal. Where is the other signal? By the way, I'm now leaving for a few hours.

[Akbar]

t*1000 means time axis in milliseconds. If i write y = wavrecord(dt*Fs, Fs,1); instead of y = wavrecord(dt*Fs, Fs, 2); i see only blue signal at the output. But when i write y = wavrecord(dt*Fs, Fs, 2); then both green and blue appear. As a reference type help wavrecord in command window. Ok, i will be here, waiting for your responce, thank you

[Akbar]

FYI when i write y = wavrecord(dt*Fs, Fs, 3); output is only empty screen with grids, because my sound card is only two channel (actually most of the soundcards on laptops are two channel ones)

[Star Strider]

If your code is:

y = wavrecord(dt*Fs, Fs, 2);

what size is y?

How are you plotting y?

[Image Analyst]

Thank you for (finally) being honest about it being one of your exam questions. What will be your solution, since you can't use Star's answer as your own on your exam? Good luck on your exam. Next time, please tag the post as homework or exam before getting solutions from us.

[Star Strider]

How interesting! I thought it was simply a home physics experiment.

It was a fun exercise for me, though it took me a couple hours to get it working correctly.

You might want to inform your professor. Had I known it was an exam question, I’d probably have written the code for fun, but not have posted it.

[Star Strider]

I’m not angry!!

I’m simply amused that your professor and university don’t care whether you learn maths, science, and programming, so long as you get an acceptable answer.

When I was an undergraduate, we had to write out and sign ‘The Pledge’ at the end of every assignment and exam. I still remember it verbatim:

• ‘I pledge on my honour that I have neither given nor received any unacknowledged aid on this [assignment, examination, etc.].’
• {Signature}

Being found in violation of ‘The Pledge’ would ‘sever your connection from the University’. Ours was — and remains — a highly-respected university, and the Honour Code was one of the reasons. We placed everyone who wanted to pursue postgraduate education in highly-respected graduate academic and professional schools. A couple went on to win Nobel Prizes, and others to other high honours.

[Akbar]

Yes, you are right I must do it myself. I just wanted to speed up the process because i'm busy with other exams, labs and projects. Star Strider, my university and professor do care about what I learn and i'm trying to learn, i'm not the one who copies h/w from internet and submits it to teacher. I really hate copying from internet without understanding anything.

Unfortunately your problem doesn't solve my problem properly. I decided to change the algorithm. I want to plot data coming from every mic. separately, then find the time of highest peak in each plot. The difference between these highest peak times will give me delay that i was looking for.

y in my code is a two column vector, each column represents data coming to each mic1 and mic2. I want to separate that data but when i say:

    y = wavrecord(dt*Fs, Fs,2) %y is two column vector here
    plot(t*1000,y);

MATLAB plots both data simultaneously as in figures that i attached in comments. when i say:

    y = wavrecord(dt*Fs, Fs,1) %y is one column vector here
    plot(t*1000,y);

It plots only one mics data, i can't get MATLAB plot the other one.

And then let's say i plot it, how do i find the peak? Can i use max() function? Or which function would do the job best?

[Star Strider]

I find the peaks of both channels in my code. I doubt if it has to be adapted to find the peaks, since it defines a peak as 95% of the maximum of the bandpass-filtered signals in each channel.

To plot both channels (assuming y is (Nx2)):

y = wavrecord(dt*Fs, Fs,2) %y is two column vector here

figure(1)
subplot(2,1,1)
plot(t*1000,y(:,1));
grid
subplot(2,1,2)
plot(t*1000,y(:,2));
grid

[Akbar]

In this code, after the end of a FOR loop Star Strider code starts. I cannot connect both pieces of code.

    clear; figure; grid on; hold on; 
Fs = 4000; % sampling rate in Hz 
dt = 4;    % duration in seconds 
% get data 
y = wavrecord(dt*Fs, Fs, 2);  % it's a vector with two columns (data from two mics)
% convert to time (sec)
tmax = length(y)/Fs; 
t = linspace(0, tmax, dt*Fs); 
% plot 
plot(t*1000,y); 
axis([500 tmax*1000 0 0.9]); 
xlabel('time (ms)'); 
% --- Analysis --- 
j = 1; 
for i=1:length(y) 
  if y(i)>0.15 
      pick(j) = 1000*i/Fs; 
      fprintf('%3d --> %8.1f ms\n',j, pick(j)); 
      j=j+1; 
  end 
end

 %Get figure and data, convert to numeric from cell:
figfile = 'Data2.fig';
ssf = openfig(figfile);
hf2 = get(gca, 'Children');
hf3xc = get(hf2, 'XData');
hf3yc = get(hf2, 'YData');
hf3xd = cell2mat(hf3xc);
hf3yd = cell2mat(hf3yc);

 % Bandpass filter:
[b,a] = butter(3, [0.2  .9])
y1f = filtfilt(b, a, hf3yd(1,:));
y2f = filtfilt(b, a, hf3yd(2,:));
pulsidx1 = find(y1f >= max(y1f)*0.5);
pulsidx2 = find(y2f >= max(y2f)*0.5);
pulse1 = pulsidx1;
pulse2 = pulsidx2;
Fs = 4000;                       % Sampling frequency (Hz)
Ts = 1/Fs;                       % Sampling interval (seconds)
% Equalise vector lengths:  
pl1 = length(pulse1);
pl2 = length(pulse2);
vlen = min(pl1, pl2);
pulsx = pulse1(1:vlen);         % Lengths in samples
pulsy = pulse2(1:vlen);
pulsxt = pulsx/Fs;              % Lengths in seconds
pulsyt = pulsy/Fs;
% Subtract:
%PulsDs = pulsx - pulsy;         % Difference in microphones arrival times (samples)
PulsDt = pulsxt - pulsyt;       % Difference in microphones arrival times (seconds)
% Take mean:
%PulsDsM = mean(abs(PulsDs));
PulsDtM = mean(abs(PulsDt));
MicDist = 1.58;                 % Distance between microphones (m)
Csc = MicDist/PulsDtM           % Calculated speed of sound (m/s)
fprintf(1,'\n\n\tSpeed of sound for %s = %.3f m/s\n\n', figfile, Csc)

[Star Strider]

You don’t need the section that gets the data from your figure files, since you already have it. You can delete or comment-out that part. I included it simply as an illustration of how to do that, and how I recovered your data.

In my code:

• The hf3xd variable are the time vectors, a (2xN) matrix of duplicate values, the same for each channel.
• The hf3yd variable are the two channels of your audio data in a (2xN) matrix the same size as hf3xd.

Assign your t and y data to my variables hf3xd and hf3yd respectively just before the ‘% Bandpass filter’ comment, and you can use it with your new data without modifying any of the rest of my code.

[Akbar]

Before Bandpass filter, where is the mistake?

    clear; figure; grid on; hold on; 
Fs = 4000; % sampling rate in Hz 
dt = 4;    % duration in seconds 
% get data 
y = wavrecord(dt*Fs, Fs, 2);  % it's a vector with two columns (data from two mics)
% convert to time (sec)
tmax = length(y)/Fs; 
t = linspace(0, tmax, dt*Fs); 
% plot 
plot(t*1000,y); 
axis([500 tmax*1000 0 0.9]); 
xlabel('time (ms)'); 
% --- Analysis --- 
j = 1; 
for i=1:length(y) 
  if y(i)>0.15 
      pick(j) = 1000*i/Fs; 
      fprintf('%3d --> %8.1f ms\n',j, pick(j)); 
      j=j+1; 
  end 
end

 hf3xd = t;
 hf3yd = y;

 % Bandpass filter:
[b,a] = butter(3, [0.2  .9])
y1f = filtfilt(b, a, hf3yd(1,:));
y2f = filtfilt(b, a, hf3yd(2,:));
pulsidx1 = find(y1f >= max(y1f)*0.5);
pulsidx2 = find(y2f >= max(y2f)*0.5);
pulse1 = pulsidx1;
pulse2 = pulsidx2;
Fs = 4000;                       % Sampling frequency (Hz)
Ts = 1/Fs;                       % Sampling interval (seconds)
 % Equalise vector lengths:  
pl1 = length(pulse1);
pl2 = length(pulse2);
vlen = min(pl1, pl2);
pulsx = pulse1(1:vlen);         % Lengths in samples
pulsy = pulse2(1:vlen);
pulsxt = pulsx/Fs;              % Lengths in seconds
pulsyt = pulsy/Fs;
% Subtract:
%PulsDs = pulsx - pulsy;         % Difference in microphones arrival times (samples)
PulsDt = pulsxt - pulsyt;       % Difference in microphones arrival times (seconds)
% Take mean:
%PulsDsM = mean(abs(PulsDs));
PulsDtM = mean(abs(PulsDt));
MicDist = 1.58;                 % Distance between microphones (m)
Csc = MicDist/PulsDtM           % Calculated speed of sound (m/s)
fprintf(1,'\n\n\tSpeed of sound for %s = %.3f m/s\n\n', figfile, Csc)

Thank you!

[Star Strider]

I can’t answer that without knowing more.

• What errors is it throwing? What lines?
• What is it doing that it shouldn’t, or what isn’t it doing that it should?
• What are the sizes of your t and y matrices?

[Akbar]

I just want to connect both pieces of code.

[Star Strider]

I think I see the problem.

Try this:

 hf3xd = t';
 hf3yd = y';

Your figures produced row vectors, but you’re collecting data as column vectors. Transpose, and you should be fine.

[Akbar]

Unfortunately it didn't help. Whatever, thank you very very much! :) Error in a last line (line52).

[Star Strider]

I assigned the name of each of your files figfile. That made it easier to load them, and then to write out the speed of sound for the appropriate figure. (See my previous posts for examples of that output.) You can either rewrite the line, or name figfile as any string variable.

One option:

figfile = datestr(now, 'yyyy-mm-dd HH.MM.SS');

That way, it prints out the date and time of the experiment, as well as the calculated speed of sound.

[triveni p]

can you just

 suggest any code for butterworth in image noise re

moval