forked from UKY-Distributed-Audio-Lab/Array-Toolbox
-
Notifications
You must be signed in to change notification settings - Fork 0
Expand file tree
/
Copy pathasa.m
More file actions
309 lines (272 loc) · 13.3 KB
/
asa.m
File metadata and controls
309 lines (272 loc) · 13.3 KB
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
function asa()
%
% This function computes an Auditory Scene Analysis (ASA) based on a
% recording from a distributed microphone system. The ASA consists of a
% text file identifying an active sound source, its location, and its
% connection to sound sources at other times. The function prompts for a
% multichannel wave file from a distributed microphone recording. It then
% prompts for a .txt file containing the microphone array information.
% It then prompts through a series of dialog boxes for pressure,
% temperature, and humidity (for sound speed calculation), and FOV regions.
% Users enter the range in the Z direction and then selects either square
% or rectangular regions using mouse clicks.
%
% It then performs Sound Source Location (SSL) using a Coherent Steered
% Response Power (CSRP) algorithm.
% CFAR threshold (1e-5) is applied to identify sound sources. A text file
% is then created to identify detected sources, strength of their detection
% statistics, and their locations in time and space. For each detection,
% a row is created with the following format:
% strms1 = [timestamp, detection statistic in excess of threshold, x, y, z postions in meters]
%
% For now this information is written as a variable in a .mat file along
% with parameters in a data structure FOVP, such as sample frequency (fs)
% and speed of sound (c).
%
% Required Functions:
% signifslist.m
%
% Written by Kevin D. Donohue ([email protected]) May 2012
% Modified by Kirstin Brangers July 2012
%
fovp.graphicout_flag = 1; % Set to 1 for collapsed acoustic images for each frame while processing
fovp.fc = 400; % High-pass cutoff freq. filter for signal conditioning
fovp.fs = 16000; % Sampling rate in Hz for resampling audio file
fovp.beta = .85; % Beta for whitening
% fovp.shape = 1.0; % Shape parameter for probability of false alarm computation
% Time window for Frequency Domain Block Processing
fovp.trez = 100e-3; % In seconds
fovp.tinc = fovp.trez/2; % Time increment in seconds
% CFAR Threshold
fovp.pcfar = 1e-2; % False alarm probability
% Room Resolution: Step through Cartesian grid for mic and sound source plane
rad = [.42, .42, .2]; % Threshold neighborhood for CFAR detection [x,y,z] in meters
rez = [.04, .04, .4]; % Grid Resolution in Meters [x,y,z] in meters
% Compute Region Extension to Compensate for Threshold Neighborhood
fovp.threshneigh = round(rad ./ rez); % Convert neighborhood in meters to samples
% Prompt for wave file:
[fn, pn] = uigetfile('*.wav','Select Multichannel Recording.');
if isnumeric(fn) % If menu cancelled, quit
error('file error')
end
% Prompt for microphone position file
[fmp, pmp] = uigetfile('*.txt','Select text file with microphone positions.');
if isnumeric(fn) % If menu cancelled, quit
error('file error')
end
fovp.mp = load([pmp, fmp]); % Load file with mic positions
% Prompt for Environment Parameters
prompt={'Temperature (degrees C):','Pressure (mmHg):','Relative Humidity (percent):'};
name='Parameters to Determine Sound Speed';
numlines=1;
defaultanswer={'24.0','29.06', '61.5'};
answer=inputdlg(prompt,name,numlines,defaultanswer,'on');
tmp = str2double(answer{1}); % Temperature
pres = str2double(answer{2}); % Pressure
rh = str2double(answer{3}); % Relative Humidity
% Compute Speed of Sound with Entered Parameters
fovp.c = SpeedOfSound(tmp,rh,pres);
% Compute Limits for Initial Plots
minx = min(fovp.mp(1,:));
maxx = max(fovp.mp(1,:));
miny = min(fovp.mp(2,:));
maxy = max(fovp.mp(2,:));
minz = min(fovp.mp(3,:));
minz = min([minz,0]);
maxz = max(fovp.mp(3,:));
% Plot Microphone Positions
figure(1)
plot3(fovp.mp(1,:),fovp.mp(2,:),fovp.mp(3,:),'xb') % Plot mic postions in blue with 'x'
xlabel('X Meters');
ylabel('Y Meters');
zlabel('Z Meters');
title(['Microphone Positions with Sound Speed = ' num2str(fovp.c) ' m/s'])
axis([minx(1)-1, maxx(1)+1, miny(1)-1, maxy(1)+1, minz(1),maxz(1)+.1])
grid on
pause(1)
% % Prompt User for Height Coordinates
% % Default coords for Z position - range from .6 meters (~2ft) to 2 meters (~6.5ft)
% % NOTE: This range is used for all regions selected in FOV
% prompt={'Enter first Z coord of FOV (meters):', 'Enter second Z coord of FOV (meters):'};
% name='Parameters to Determine Height of FOV Region';
% numlines=1;
% defaultanswer={'0.6','2.0'}; % Default Z range
% answer=inputdlg(prompt,name,numlines,defaultanswer,'on');
% z1 = str2double(answer{1}); % First coord for Z range
% z2 = str2double(answer{2}); % Second coord for Z range
% 2D Plot So X,Y Coordinates Can Be Selected with Mouse Click
figure(1)
plot(fovp.mp(1,:),fovp.mp(2,:),'xb') % Plot mic postions in blue with 'x'
xlabel('X Meters');
ylabel('Y Meters');
title(['Microphone Positions with Sound Speed = ' num2str(fovp.c) ' m/s'])
axis([minx(1)-1, maxx(1)+1, miny(1)-1, maxy(1)+1])
grid on
pause(1)
% Prompt User for Region(s) of Interest
fovflag = 'Yes'; % Turn on Fov Flag to prompt user for initial region
fovcnt = 0; % Initialize counter for number of FOV lines
xyposS = 1; % Initialize counter to hold X,Y coordinates for selecting source
while fovflag(1) == 'Y'
% Prompt user to either select a source with radius or region of interest
roi = questdlg('Select Source or Region.', 'Region of Interest', 'Source' ,'Region', 'Source');
% Prompt User for Height Coordinates
% Default coords for Z position - range from .6 meters (~2ft) to 2 meters (~6.5ft)
prompt={'Enter first Z coord of FOV (meters):', 'Enter second Z coord of FOV (meters):'};
name='Parameters to Determine Height of FOV Region';
numlines=1;
defaultanswer={'0.6','2.0'}; % Default Z range
answer=inputdlg(prompt,name,numlines,defaultanswer,'on');
z1 = str2double(answer{1}); % First coord for Z range
z2 = str2double(answer{2}); % Second coord for Z range
fovm(fovcnt+1,3)=z1;
fovm(fovcnt+2,3)=z2;
if roi(1) == 'S' % If user selects 'Source'
% Select Source of Interest - One Mouse Click
h = helpdlg('Select Position of Sound Source with Crosshairs','Select Source Position.');
uiwait(h); % Delete dialog box when user selects OK
[xCoord(xyposS),yCoord(xyposS)] = ginput(1); % Get X,Y position of mouse click
% Ask user for radius around source
rdef = 0.0; % Default radius
prompt={'Enter radius around source (meters):'};
name='Parameters to Determine Radius of FOV Region';
numlines=1;
defaultanswer={num2str(rdef,3)};
answer=inputdlg(prompt,name,numlines,defaultanswer);
radxy = str2double(answer{1}); % Radius entered by user
% If radius is less than default radius, change to default radius
if radxy < rdef
radxy = rdef;
end
% Radius cannot be negative
if radxy < 0
radxy = 0;
end
% Add radius to selected source position to obtain opposite corners
fovm(fovcnt+1,1)= xCoord(xyposS)-radxy; % Corner 1
fovm(fovcnt+1,2)= yCoord(xyposS)-radxy;
fovm(fovcnt+2,1)= xCoord(xyposS)+radxy; % Corner 2
fovm(fovcnt+2,2)= yCoord(xyposS)+radxy;
% Get values for graphing purposes
x = min(fovm(fovcnt+1,1),fovm(fovcnt+2,1)); % X position
y = min(fovm(fovcnt+1,2),fovm(fovcnt+2,2)); % Y position
w = abs((fovm(fovcnt+1,1)-fovm(fovcnt+2,1))); % Width of rect
h = abs((fovm(fovcnt+1,2)-fovm(fovcnt+2,2))); % Height of rect
% Plot selected source including radius
figure(1)
hold on
h = rectangle('Position',[x, y, w, h],'EdgeColor','r');
hold off
else % If user selects 'Region'
% Select Region of Interest - Two Mouse Clicks
h = helpdlg('Select Region of Interest with Crosshairs (opposite corners)','Select Source Corners.');
uiwait(h);
[xCoordA,yCoordA] = ginput(1); % First selected corner
% Show first selected corner on plot
figure(1)
hold on
h = plot(xCoordA,yCoordA,'r.-'); % Plot selected corner in red with '.'
hold off
% Get second corner
[xCoordB,yCoordB] = ginput(1);
% Opposite corners
fovm(fovcnt+1,1)= xCoordA; % Corner 1 X1
fovm(fovcnt+1,2)= yCoordA; % Y1
fovm(fovcnt+2,1)= xCoordB; % Corner 2 X2
fovm(fovcnt+2,2)= yCoordB; % Y2
% Get smaller value for graphing purposes
x = min(xCoordA,xCoordB); % X position
y = min(yCoordA,yCoordB); % Y position
w = abs(xCoordA-xCoordB); % Width of rectangle
h = abs(yCoordA-yCoordB); % Height of rectangle
% Plot selected region
figure(1)
hold on
h = rectangle('Position',[x, y, w, h],'EdgeColor','r');
hold off
end
% Increment Counters
xyposS = xyposS+1;
fovcnt = fovcnt+2;
fovflag = questdlg('Add another block of interest?', 'Field of View');
if fovflag(1) == 'C' % If cancel is selected, quit
error('User exited program.')
end
end
% Delete previously chosen regions on 2d plot
delete(h);
% Display Selected Corners in Matlab Command Window
lenfovcnt = fovcnt;
fovcnt = 0;
region = 1;
xyposS = 1;
fprintf('The selected regions of interest: \n');
while lenfovcnt ~= 0
fprintf(['The Source of Interest for Region ',num2str(region), ' is at corners: ( ',num2str(fovm(fovcnt+1,:)), ' ) and ( ' ,num2str(fovm(fovcnt+2,:)),' ) \n \n']);
lenfovcnt = lenfovcnt-2;
fovcnt = fovcnt+2;
region = region+1;
xyposS = xyposS+1;
end
% Use FOV points to create grid for analysis
% Plot user selected points only - no CFAR Threshold
for k=1:2:fovcnt
for kc=1:3
xmm = min(fovm(k:k+1,kc));
xmx = max(fovm(k:k+1,kc));
fov(kc,1) = xmm(1);
fov(kc,2) = xmx(1);
end
vertx{fix((k-1)/2+1)} = [fov(:,1), [fov(1:2,1); fov(3,2)], [fov(1,1); fov(2,2); fov(3,1)], [fov(1,1); fov(2:3,2)], ...
[fov(1,2); fov(2:3,1)], [fov(1,2); fov(2,1); fov(3,2)], [fov(1:2,2); fov(3,1)], fov(:,2)];
fovp.sgrid{fix((k-1)/2+1)} = {[fov(1,1):rez(1):fov(1,2)], [fov(2,1):rez(2):fov(2,2)], [fov(3,1):rez(3):fov(3,2)]};
end
fnum = length(vertx); % Number of FOV region blocks
% Plot mics with FOV regions
figure(1)
plot3(fovp.mp(1,:),fovp.mp(2,:),fovp.mp(3,:),'xb')
xlabel('X Meters');
ylabel('Y Meters');
zlabel('Z Meters');
title(['Microphone Positions with Sound Speed = ' num2str(fovp.c) ' m/s'])
grid on
pause(1)
hold on
% Skech layout of microphone and FOV regions before adding CFAR Threshold
for kv=1:fnum
% Y lines
plot3(vertx{kv}(1,1:2),vertx{kv}(2,1:2),vertx{kv}(3,1:2),'k.-')
plot3(vertx{kv}(1,3:4),vertx{kv}(2,3:4),vertx{kv}(3,3:4),'k.-')
plot3(vertx{kv}(1,5:6),vertx{kv}(2,5:6),vertx{kv}(3,5:6),'k.-')
plot3(vertx{kv}(1,7:8),vertx{kv}(2,7:8),vertx{kv}(3,7:8),'k.-')
% Z lines
plot3(vertx{kv}(1,1:2:3),vertx{kv}(2,1:2:3),vertx{kv}(3,1:2:3),'k.-')
plot3(vertx{kv}(1,2:2:4),vertx{kv}(2,2:2:4),vertx{kv}(3,2:2:4),'k.-')
plot3(vertx{kv}(1,5:2:7),vertx{kv}(2,5:2:7),vertx{kv}(3,5:2:7),'k.-')
plot3(vertx{kv}(1,6:2:8),vertx{kv}(2,6:2:8),vertx{kv}(3,6:2:8),'k.-')
% X lines
plot3(vertx{kv}(1,1:4:5),vertx{kv}(2,1:4:5),vertx{kv}(3,1:4:5),'k.-')
plot3(vertx{kv}(1,2:4:6),vertx{kv}(2,2:4:6),vertx{kv}(3,2:4:6),'k.-')
plot3(vertx{kv}(1,3:4:7),vertx{kv}(2,3:4:7),vertx{kv}(3,3:4:7),'k.-')
plot3(vertx{kv}(1,4:4:8),vertx{kv}(2,4:4:8),vertx{kv}(3,4:4:8),'k.-')
end
hold off
pause(.1) % Pause to let graphics operate
% Add on CFAR Threshold Detection
for k=1:2:fovcnt
for kc=1:3
xmm = min(fovm(k:k+1,kc));
xmx = max(fovm(k:k+1,kc));
fov(kc,1) = xmm(1)-fovp.threshneigh(kc)*rez(kc); % Add extra to account for CFAR averaging
fov(kc,2) = xmx(1)+fovp.threshneigh(kc)*rez(kc); % Add extra to account for CFAR averaging
end
vertx{fix((k-1)/2+1)} = [fov(:,1), [fov(1:2,1); fov(3,2)], [fov(1,1); fov(2,2); fov(3,1)], [fov(1,1); fov(2:3,2)], ...
[fov(1,2); fov(2:3,1)], [fov(1,2); fov(2,1); fov(3,2)], [fov(1:2,2); fov(3,1)], fov(:,2)];
fovp.sgrid{fix((k-1)/2+1)} = {[fov(1,1):rez(1):fov(1,2)], [fov(2,1):rez(2):fov(2,2)], [fov(3,1):rez(3):fov(3,2)]};
end
% Assign wavefile recording filename to input data structure
fovp.fn = fn; % file name
fovp.pn = pn; % path make empty [] if in current working directory
% Detection significant sound sources in each frame
strms1 = signifslist(fovp); % Find sources in selected regions
save([fn(1:end-4) '.mat'],'strms1', 'fovp') % Save scipt of all detected sources