[READ-ONLY] Mirror of https://github.com/jmrplens/EASEPostFile2Matlab. Representa los datos de historia temporal obtenidos con EASE 4 y calcula la teoría revisada corregida.
Add files via upload
authored by
Jose Manuel Requena Plens
and committed by
GitHub
(Jun 3, 2018, 2:58 AM +0200)
d8e3f57b
cc483730
+3507
+58
include/EASEPost2Distan.m
+58
include/EASEPost2Distan.m
···
1
+
function Centro01 = EASEPost2Distan(filename, startRow, endRow)
2
+
%IMPORTFILE Import numeric data from a text file as a matrix.
3
+
% CENTRO01 = IMPORTFILE(FILENAME) Reads data from text file FILENAME for
4
+
% the default selection.
5
+
%
6
+
% CENTRO01 = IMPORTFILE(FILENAME, STARTROW, ENDROW) Reads data from rows
7
+
% STARTROW through ENDROW of text file FILENAME.
8
+
%
9
+
% Example:
10
+
% Centro01 = importfile('Centro_01.txt', [16,22], [16,22]);
11
+
%
12
+
% See also TEXTSCAN.
13
+
14
+
% Auto-generated by MATLAB on 2018/02/28 00:31:51
15
+
16
+
%% Initialize variables.
17
+
delimiter = ' ';
18
+
if nargin<=2
19
+
startRow = [16,22];
20
+
endRow = [16,22];
21
+
end
22
+
23
+
%% Format for each line of text:
24
+
% column3: double (%f)
25
+
% column4: double (%f)
26
+
% column5: double (%f)
27
+
% For more information, see the TEXTSCAN documentation.
28
+
formatSpec = '%*s%*s%f%f%f%*s%*s%*s%*s%*s%*s%*s%*s%*s%*s%*s%*s%*s%*s%*s%*s%*s%*s%*s%*s%*s%[^\n\r]';
29
+
30
+
%% Open the text file.
31
+
fileID = fopen(filename,'r');
32
+
33
+
%% Read columns of data according to the format.
34
+
% This call is based on the structure of the file used to generate this
35
+
% code. If an error occurs for a different file, try regenerating the code
36
+
% from the Import Tool.
37
+
textscan(fileID, '%[^\n\r]', startRow(1)-1, 'WhiteSpace', '', 'ReturnOnError', false);
38
+
dataArray = textscan(fileID, formatSpec, endRow(1)-startRow(1)+1, 'Delimiter', delimiter, 'MultipleDelimsAsOne', true, 'TextType', 'string', 'ReturnOnError', false, 'EndOfLine', '\r\n');
39
+
for block=2:length(startRow)
40
+
frewind(fileID);
41
+
textscan(fileID, '%[^\n\r]', startRow(block)-1, 'WhiteSpace', '', 'ReturnOnError', false);
42
+
dataArrayBlock = textscan(fileID, formatSpec, endRow(block)-startRow(block)+1, 'Delimiter', delimiter, 'MultipleDelimsAsOne', true, 'TextType', 'string', 'ReturnOnError', false, 'EndOfLine', '\r\n');
43
+
for col=1:length(dataArray)
44
+
dataArray{col} = [dataArray{col};dataArrayBlock{col}];
45
+
end
46
+
end
47
+
48
+
%% Close the text file.
49
+
fclose(fileID);
50
+
51
+
%% Post processing for unimportable data.
52
+
% No unimportable data rules were applied during the import, so no post
53
+
% processing code is included. To generate code which works for
54
+
% unimportable data, select unimportable cells in a file and regenerate the
55
+
% script.
56
+
57
+
%% Create output variable
58
+
Centro01 = [dataArray{1:end-1}];
+55
include/EASEPost2Info.m
+55
include/EASEPost2Info.m
···
1
+
function Centro01 = EASEPost2Info(filename, startRow, endRow)
2
+
%IMPORTFILE Import numeric data from a text file as a matrix.
3
+
% CENTRO01 = IMPORTFILE(FILENAME) Reads data from text file FILENAME for
4
+
% the default selection.
5
+
%
6
+
% CENTRO01 = IMPORTFILE(FILENAME, STARTROW, ENDROW) Reads data from rows
7
+
% STARTROW through ENDROW of text file FILENAME.
8
+
%
9
+
% Example:
10
+
% Centro01 = importfile('Centro_01.txt', [15,21], [15,21]);
11
+
%
12
+
% See also TEXTSCAN.
13
+
14
+
% Auto-generated by MATLAB on 2018/02/28 00:51:16
15
+
16
+
%% Initialize variables.
17
+
delimiter = ' ';
18
+
if nargin<=2
19
+
startRow = [15,21];
20
+
endRow = [15,21];
21
+
end
22
+
23
+
%% Format for each line of text:
24
+
% column2: text (%s)
25
+
% For more information, see the TEXTSCAN documentation.
26
+
formatSpec = '%*s%s%*s%*s%*s%*s%*s%*s%*s%*s%*s%*s%*s%*s%*s%*s%*s%*s%*s%*s%*s%*s%*s%*s%*s%*s%[^\n\r]';
27
+
28
+
%% Open the text file.
29
+
fileID = fopen(filename,'r');
30
+
31
+
%% Read columns of data according to the format.
32
+
% This call is based on the structure of the file used to generate this
33
+
% code. If an error occurs for a different file, try regenerating the code
34
+
% from the Import Tool.
35
+
textscan(fileID, '%[^\n\r]', startRow(1)-1, 'WhiteSpace', '', 'ReturnOnError', false);
36
+
dataArray = textscan(fileID, formatSpec, endRow(1)-startRow(1)+1, 'Delimiter', delimiter, 'MultipleDelimsAsOne', true, 'TextType', 'string', 'ReturnOnError', false, 'EndOfLine', '\r\n');
37
+
for block=2:length(startRow)
38
+
frewind(fileID);
39
+
textscan(fileID, '%[^\n\r]', startRow(block)-1, 'WhiteSpace', '', 'ReturnOnError', false);
40
+
dataArrayBlock = textscan(fileID, formatSpec, endRow(block)-startRow(block)+1, 'Delimiter', delimiter, 'MultipleDelimsAsOne', true, 'TextType', 'string', 'ReturnOnError', false, 'EndOfLine', '\r\n');
41
+
dataArray{1} = [dataArray{1};dataArrayBlock{1}];
42
+
end
43
+
44
+
%% Close the text file.
45
+
fclose(fileID);
46
+
47
+
%% Post processing for unimportable data.
48
+
% No unimportable data rules were applied during the import, so no post
49
+
% processing code is included. To generate code which works for
50
+
% unimportable data, select unimportable cells in a file and regenerate the
51
+
% script.
52
+
53
+
%% Create output variable
54
+
Centro01 = [dataArray{1:end-1}];
55
+
+77
include/EASEPost2Matrix.m
+77
include/EASEPost2Matrix.m
···
1
+
function Centro1 = EASEPost2Matrix(filename, startRow, endRow)
2
+
%IMPORTFILE Import numeric data from a text file as a matrix.
3
+
% CENTRO1 = IMPORTFILE(FILENAME) Reads data from text file FILENAME for
4
+
% the default selection.
5
+
%
6
+
% CENTRO1 = IMPORTFILE(FILENAME, STARTROW, ENDROW) Reads data from rows
7
+
% STARTROW through ENDROW of text file FILENAME.
8
+
%
9
+
% Example:
10
+
% Centro1 = importfile('Centro_1.txt', 31, 61966);
11
+
%
12
+
% See also TEXTSCAN.
13
+
14
+
% Auto-generated by MATLAB on 2018/02/27 22:06:29
15
+
16
+
%% Initialize variables.
17
+
delimiter = ' ';
18
+
if nargin<=2
19
+
startRow = 31;
20
+
endRow = inf;
21
+
end
22
+
23
+
%% Format for each line of text:
24
+
% column2: double (%f)
25
+
% column3: double (%f)
26
+
% column4: double (%f)
27
+
% column5: double (%f)
28
+
% column6: double (%f)
29
+
% column7: double (%f)
30
+
% column8: double (%f)
31
+
% column9: double (%f)
32
+
% column10: double (%f)
33
+
% column11: double (%f)
34
+
% column12: double (%f)
35
+
% column13: double (%f)
36
+
% column14: double (%f)
37
+
% column15: double (%f)
38
+
% column16: double (%f)
39
+
% column17: double (%f)
40
+
% column18: double (%f)
41
+
% column19: double (%f)
42
+
% column20: double (%f)
43
+
% column21: double (%f)
44
+
% column22: double (%f)
45
+
% column23: double (%f)
46
+
% For more information, see the TEXTSCAN documentation.
47
+
formatSpec = '%*s%f%f%f%f%f%f%f%f%f%f%f%f%f%f%f%f%f%f%f%f%f%f%*s%*s%*s%[^\n\r]';
48
+
49
+
%% Open the text file.
50
+
fileID = fopen(filename,'r');
51
+
52
+
%% Read columns of data according to the format.
53
+
% This call is based on the structure of the file used to generate this
54
+
% code. If an error occurs for a different file, try regenerating the code
55
+
% from the Import Tool.
56
+
textscan(fileID, '%[^\n\r]', startRow(1)-1, 'WhiteSpace', '', 'ReturnOnError', false);
57
+
dataArray = textscan(fileID, formatSpec, endRow(1)-startRow(1)+1, 'Delimiter', delimiter, 'MultipleDelimsAsOne', true, 'TextType', 'string', 'EmptyValue', NaN, 'ReturnOnError', false, 'EndOfLine', '\r\n');
58
+
for block=2:length(startRow)
59
+
frewind(fileID);
60
+
textscan(fileID, '%[^\n\r]', startRow(block)-1, 'WhiteSpace', '', 'ReturnOnError', false);
61
+
dataArrayBlock = textscan(fileID, formatSpec, endRow(block)-startRow(block)+1, 'Delimiter', delimiter, 'MultipleDelimsAsOne', true, 'TextType', 'string', 'EmptyValue', NaN, 'ReturnOnError', false, 'EndOfLine', '\r\n');
62
+
for col=1:length(dataArray)
63
+
dataArray{col} = [dataArray{col};dataArrayBlock{col}];
64
+
end
65
+
end
66
+
67
+
%% Close the text file.
68
+
fclose(fileID);
69
+
70
+
%% Post processing for unimportable data.
71
+
% No unimportable data rules were applied during the import, so no post
72
+
% processing code is included. To generate code which works for
73
+
% unimportable data, select unimportable cells in a file and regenerate the
74
+
% script.
75
+
76
+
%% Create output variable
77
+
Centro1 = [dataArray{1:end-1}];
+54
include/EASEPost2Otros.m
+54
include/EASEPost2Otros.m
···
1
+
function Esquina1 = EASEPost2Otros(filename, startRow, endRow)
2
+
%IMPORTFILE Import numeric data from a text file as a matrix.
3
+
% ESQUINA1 = IMPORTFILE(FILENAME) Reads data from text file FILENAME for
4
+
% the default selection.
5
+
%
6
+
% ESQUINA1 = IMPORTFILE(FILENAME, STARTROW, ENDROW) Reads data from rows
7
+
% STARTROW through ENDROW of text file FILENAME.
8
+
%
9
+
% Example:
10
+
% Esquina1 = importfile('Esquina_1.txt', [18,24,29], [18,24,29]);
11
+
%
12
+
% See also TEXTSCAN.
13
+
14
+
% Auto-generated by MATLAB on 2018/03/02 01:50:48
15
+
16
+
%% Initialize variables.
17
+
delimiter = ' ';
18
+
if nargin<=2
19
+
startRow = [18,24,29];
20
+
endRow = [18,24,29];
21
+
end
22
+
23
+
%% Format for each line of text:
24
+
% column4: double (%f)
25
+
% For more information, see the TEXTSCAN documentation.
26
+
formatSpec = '%*s%*s%*s%f%*s%*s%*s%*s%*s%*s%*s%*s%*s%*s%*s%*s%*s%*s%*s%*s%*s%*s%*s%*s%*s%*s%[^\n\r]';
27
+
28
+
%% Open the text file.
29
+
fileID = fopen(filename,'r');
30
+
31
+
%% Read columns of data according to the format.
32
+
% This call is based on the structure of the file used to generate this
33
+
% code. If an error occurs for a different file, try regenerating the code
34
+
% from the Import Tool.
35
+
textscan(fileID, '%[^\n\r]', startRow(1)-1, 'WhiteSpace', '', 'ReturnOnError', false);
36
+
dataArray = textscan(fileID, formatSpec, endRow(1)-startRow(1)+1, 'Delimiter', delimiter, 'MultipleDelimsAsOne', true, 'TextType', 'string', 'ReturnOnError', false, 'EndOfLine', '\r\n');
37
+
for block=2:length(startRow)
38
+
frewind(fileID);
39
+
textscan(fileID, '%[^\n\r]', startRow(block)-1, 'WhiteSpace', '', 'ReturnOnError', false);
40
+
dataArrayBlock = textscan(fileID, formatSpec, endRow(block)-startRow(block)+1, 'Delimiter', delimiter, 'MultipleDelimsAsOne', true, 'TextType', 'string', 'ReturnOnError', false, 'EndOfLine', '\r\n');
41
+
dataArray{1} = [dataArray{1};dataArrayBlock{1}];
42
+
end
43
+
44
+
%% Close the text file.
45
+
fclose(fileID);
46
+
47
+
%% Post processing for unimportable data.
48
+
% No unimportable data rules were applied during the import, so no post
49
+
% processing code is included. To generate code which works for
50
+
% unimportable data, select unimportable cells in a file and regenerate the
51
+
% script.
52
+
53
+
%% Create output variable
54
+
Esquina1 = [dataArray{1:end-1}];
+76
include/EASEPost2Param.m
+76
include/EASEPost2Param.m
···
1
+
function Esquina1 = EASEPost2Param(filename, startRow, endRow)
2
+
%IMPORTFILE Import numeric data from a text file as a matrix.
3
+
% ESQUINA1 = IMPORTFILE(FILENAME) Reads data from text file FILENAME for
4
+
% the default selection.
5
+
%
6
+
% ESQUINA1 = IMPORTFILE(FILENAME, STARTROW, ENDROW) Reads data from rows
7
+
% STARTROW through ENDROW of text file FILENAME.
8
+
%
9
+
% Example:
10
+
% Esquina1 = importfile('Esquina_1.txt', [11,28], [12,28]);
11
+
%
12
+
% See also TEXTSCAN.
13
+
14
+
% Auto-generated by MATLAB on 2018/03/02 02:08:49
15
+
16
+
%% Initialize variables.
17
+
delimiter = ' ';
18
+
if nargin<=2
19
+
startRow = [11,28];
20
+
endRow = [12,28];
21
+
end
22
+
23
+
%% Format for each line of text:
24
+
% column3: double (%f)
25
+
% column4: double (%f)
26
+
% column5: double (%f)
27
+
% column6: double (%f)
28
+
% column7: double (%f)
29
+
% column8: double (%f)
30
+
% column9: double (%f)
31
+
% column10: double (%f)
32
+
% column11: double (%f)
33
+
% column12: double (%f)
34
+
% column13: double (%f)
35
+
% column14: double (%f)
36
+
% column15: double (%f)
37
+
% column16: double (%f)
38
+
% column17: double (%f)
39
+
% column18: double (%f)
40
+
% column19: double (%f)
41
+
% column20: double (%f)
42
+
% column21: double (%f)
43
+
% column22: double (%f)
44
+
% column23: double (%f)
45
+
% For more information, see the TEXTSCAN documentation.
46
+
formatSpec = '%*s%*s%f%f%f%f%f%f%f%f%f%f%f%f%f%f%f%f%f%f%f%f%f%*s%*s%*s%[^\n\r]';
47
+
48
+
%% Open the text file.
49
+
fileID = fopen(filename,'r');
50
+
51
+
%% Read columns of data according to the format.
52
+
% This call is based on the structure of the file used to generate this
53
+
% code. If an error occurs for a different file, try regenerating the code
54
+
% from the Import Tool.
55
+
textscan(fileID, '%[^\n\r]', startRow(1)-1, 'WhiteSpace', '', 'ReturnOnError', false);
56
+
dataArray = textscan(fileID, formatSpec, endRow(1)-startRow(1)+1, 'Delimiter', delimiter, 'MultipleDelimsAsOne', true, 'TextType', 'string', 'ReturnOnError', false, 'EndOfLine', '\r\n');
57
+
for block=2:length(startRow)
58
+
frewind(fileID);
59
+
textscan(fileID, '%[^\n\r]', startRow(block)-1, 'WhiteSpace', '', 'ReturnOnError', false);
60
+
dataArrayBlock = textscan(fileID, formatSpec, endRow(block)-startRow(block)+1, 'Delimiter', delimiter, 'MultipleDelimsAsOne', true, 'TextType', 'string', 'ReturnOnError', false, 'EndOfLine', '\r\n');
61
+
for col=1:length(dataArray)
62
+
dataArray{col} = [dataArray{col};dataArrayBlock{col}];
63
+
end
64
+
end
65
+
66
+
%% Close the text file.
67
+
fclose(fileID);
68
+
69
+
%% Post processing for unimportable data.
70
+
% No unimportable data rules were applied during the import, so no post
71
+
% processing code is included. To generate code which works for
72
+
% unimportable data, select unimportable cells in a file and regenerate the
73
+
% script.
74
+
75
+
%% Create output variable
76
+
Esquina1 = [dataArray{1:end-1}];
+76
include/EASEPost2Param2.m
+76
include/EASEPost2Param2.m
···
1
+
function Esquina1 = EASEPost2Param2(filename, startRow, endRow)
2
+
%IMPORTFILE Import numeric data from a text file as a matrix.
3
+
% ESQUINA1 = IMPORTFILE(FILENAME) Reads data from text file FILENAME for
4
+
% the default selection.
5
+
%
6
+
% ESQUINA1 = IMPORTFILE(FILENAME, STARTROW, ENDROW) Reads data from rows
7
+
% STARTROW through ENDROW of text file FILENAME.
8
+
%
9
+
% Example:
10
+
% Esquina1 = importfile('Esquina_1.txt', [13,27], [13,27]);
11
+
%
12
+
% See also TEXTSCAN.
13
+
14
+
% Auto-generated by MATLAB on 2018/03/02 02:11:55
15
+
16
+
%% Initialize variables.
17
+
delimiter = ' ';
18
+
if nargin<=2
19
+
startRow = [13,27];
20
+
endRow = [13,27];
21
+
end
22
+
23
+
%% Format for each line of text:
24
+
% column4: double (%f)
25
+
% column5: double (%f)
26
+
% column6: double (%f)
27
+
% column7: double (%f)
28
+
% column8: double (%f)
29
+
% column9: double (%f)
30
+
% column10: double (%f)
31
+
% column11: double (%f)
32
+
% column12: double (%f)
33
+
% column13: double (%f)
34
+
% column14: double (%f)
35
+
% column15: double (%f)
36
+
% column16: double (%f)
37
+
% column17: double (%f)
38
+
% column18: double (%f)
39
+
% column19: double (%f)
40
+
% column20: double (%f)
41
+
% column21: double (%f)
42
+
% column22: double (%f)
43
+
% column23: double (%f)
44
+
% column24: double (%f)
45
+
% For more information, see the TEXTSCAN documentation.
46
+
formatSpec = '%*s%*s%*s%f%f%f%f%f%f%f%f%f%f%f%f%f%f%f%f%f%f%f%f%f%*s%*s%[^\n\r]';
47
+
48
+
%% Open the text file.
49
+
fileID = fopen(filename,'r');
50
+
51
+
%% Read columns of data according to the format.
52
+
% This call is based on the structure of the file used to generate this
53
+
% code. If an error occurs for a different file, try regenerating the code
54
+
% from the Import Tool.
55
+
textscan(fileID, '%[^\n\r]', startRow(1)-1, 'WhiteSpace', '', 'ReturnOnError', false);
56
+
dataArray = textscan(fileID, formatSpec, endRow(1)-startRow(1)+1, 'Delimiter', delimiter, 'MultipleDelimsAsOne', true, 'TextType', 'string', 'ReturnOnError', false, 'EndOfLine', '\r\n');
57
+
for block=2:length(startRow)
58
+
frewind(fileID);
59
+
textscan(fileID, '%[^\n\r]', startRow(block)-1, 'WhiteSpace', '', 'ReturnOnError', false);
60
+
dataArrayBlock = textscan(fileID, formatSpec, endRow(block)-startRow(block)+1, 'Delimiter', delimiter, 'MultipleDelimsAsOne', true, 'TextType', 'string', 'ReturnOnError', false, 'EndOfLine', '\r\n');
61
+
for col=1:length(dataArray)
62
+
dataArray{col} = [dataArray{col};dataArrayBlock{col}];
63
+
end
64
+
end
65
+
66
+
%% Close the text file.
67
+
fclose(fileID);
68
+
69
+
%% Post processing for unimportable data.
70
+
% No unimportable data rules were applied during the import, so no post
71
+
% processing code is included. To generate code which works for
72
+
% unimportable data, select unimportable cells in a file and regenerate the
73
+
% script.
74
+
75
+
%% Create output variable
76
+
Esquina1 = [dataArray{1:end-1}];
+54
include/EASEPost2Sala.m
+54
include/EASEPost2Sala.m
···
1
+
function Esquina1 = EASEPost2Sala(filename, startRow, endRow)
2
+
%IMPORTFILE Import numeric data from a text file as a matrix.
3
+
% ESQUINA1 = IMPORTFILE(FILENAME) Reads data from text file FILENAME for
4
+
% the default selection.
5
+
%
6
+
% ESQUINA1 = IMPORTFILE(FILENAME, STARTROW, ENDROW) Reads data from rows
7
+
% STARTROW through ENDROW of text file FILENAME.
8
+
%
9
+
% Example:
10
+
% Esquina1 = importfile('Esquina_1.txt', 4, 7);
11
+
%
12
+
% See also TEXTSCAN.
13
+
14
+
% Auto-generated by MATLAB on 2018/03/02 01:05:38
15
+
16
+
%% Initialize variables.
17
+
delimiter = ' ';
18
+
if nargin<=2
19
+
startRow = 4;
20
+
endRow = 7;
21
+
end
22
+
23
+
%% Format for each line of text:
24
+
% column2: double (%f)
25
+
% For more information, see the TEXTSCAN documentation.
26
+
formatSpec = '%*s%f%*s%*s%*s%*s%*s%*s%*s%*s%*s%*s%*s%*s%*s%*s%*s%*s%*s%*s%*s%*s%*s%*s%*s%*s%[^\n\r]';
27
+
28
+
%% Open the text file.
29
+
fileID = fopen(filename,'r');
30
+
31
+
%% Read columns of data according to the format.
32
+
% This call is based on the structure of the file used to generate this
33
+
% code. If an error occurs for a different file, try regenerating the code
34
+
% from the Import Tool.
35
+
textscan(fileID, '%[^\n\r]', startRow(1)-1, 'WhiteSpace', '', 'ReturnOnError', false);
36
+
dataArray = textscan(fileID, formatSpec, endRow(1)-startRow(1)+1, 'Delimiter', delimiter, 'MultipleDelimsAsOne', true, 'TextType', 'string', 'ReturnOnError', false, 'EndOfLine', '\r\n');
37
+
for block=2:length(startRow)
38
+
frewind(fileID);
39
+
textscan(fileID, '%[^\n\r]', startRow(block)-1, 'WhiteSpace', '', 'ReturnOnError', false);
40
+
dataArrayBlock = textscan(fileID, formatSpec, endRow(block)-startRow(block)+1, 'Delimiter', delimiter, 'MultipleDelimsAsOne', true, 'TextType', 'string', 'ReturnOnError', false, 'EndOfLine', '\r\n');
41
+
dataArray{1} = [dataArray{1};dataArrayBlock{1}];
42
+
end
43
+
44
+
%% Close the text file.
45
+
fclose(fileID);
46
+
47
+
%% Post processing for unimportable data.
48
+
% No unimportable data rules were applied during the import, so no post
49
+
% processing code is included. To generate code which works for
50
+
% unimportable data, select unimportable cells in a file and regenerate the
51
+
% script.
52
+
53
+
%% Create output variable
54
+
Esquina1 = [dataArray{1:end-1}];
+26
include/guardarErrores.m
+26
include/guardarErrores.m
···
1
+
2
+
function guardarErrores(Ucoef,Pcoef,Distancia,Dist_P,Dist_U)
3
+
4
+
a = confint(Ucoef);
5
+
b = confint(Pcoef);
6
+
7
+
% Resolucion de errores
8
+
x = linspace(min(Distancia),max(Distancia),8);
9
+
10
+
% Puntos de error del campo util
11
+
U = Ucoef.a*x.^(Ucoef.b);
12
+
Utop = a(2)*x.^(a(4))-U;
13
+
Udow = U-a(1)*x.^(a(3));
14
+
TUerrores = table(x',U',Utop',Udow','VariableNames',{'x','y','ymas','ymenos'});
15
+
16
+
writetable(TUerrores,'UtilErrores.dat','Delimiter','\t')
17
+
18
+
% Puntos campo perjudicial
19
+
% Puntos de error del campo util
20
+
P = Pcoef.p1*x+Pcoef.p2;
21
+
Ptop = (b(2)*x+(b(4)))-P;
22
+
Pdow = P-(b(1)*x+(b(3)));
23
+
TPerrores = table(x',P',Ptop',Pdow','VariableNames',{'x','y','ymas','ymenos'});
24
+
25
+
writetable(TPerrores,'PerjudicialErrores.dat','Delimiter','\t')
26
+
end
+208
include/importarPostFiles.m
+208
include/importarPostFiles.m
···
1
+
function importarPostFiles(hObject,handles)
2
+
% Almacena en la aplicacion (handles.SPLm) la variables
3
+
% Es una variable tipo cell que en cada fila almacena los datos extraidos
4
+
% de cada fila. Si hay 20 archivos, el cell 'SPLm' tendr� 20 filas.
5
+
% En las ultimas lineas se puede observar que contiene cada columna.
6
+
%
7
+
% Distancia -> Vector de distancias entre cada receptor y la fuente
8
+
% PosReceptor -> Vector con la posicion del receptor (X,Y,Z)
9
+
% PosFuente -> Vector con la posicion de la fuente (X,Y,Z)
10
+
% Fuente -> Vector con la ID de fuente recibida cada receptor
11
+
% Receptor -> Vector con las ID de cada receptor
12
+
% Volumen -> Variable que contiene el valor del volumen de la sala
13
+
% Superficie -> Variable que contiene el valor de la superficie
14
+
% Temperatura -> Variable que contiene el valor de la temperatura
15
+
% Humedad -> Variable que contiene el valor de la humedad
16
+
% SPLm -> Cell que almacena la historia temporal de cada receptor
17
+
% RayoDirecto -> Tiempo de llegada del primer rayo a cada receptor
18
+
% RetFuente -> Retardo de la fuente
19
+
% NumRayos -> Numero de rayos recibidos por cada receptor
20
+
% Reverberacion -> Tiempo de reverberacion de cada archivo
21
+
% Absorcion -> Coeficiente de absorcion de cada archivo
22
+
% Directividad -> Directividad de las fuentes
23
+
% RuidoFondo -> Ruido de fondo de cada archivo
24
+
% SPLFuente -> Nivel SPL de las fuentes
25
+
26
+
%% Importacion de datos
27
+
ListaArchivos = handles.ListaArchivos;
28
+
29
+
%multiWaitbar(handles,'Importando archivos',N/size(PARAMtxt,1),'color','b');
30
+
% Inicializacion de variables
31
+
Distancia = zeros(size(ListaArchivos,2),1);
32
+
PosReceptor = zeros(size(ListaArchivos,2),3);
33
+
PosFuente = zeros(size(ListaArchivos,2),3);
34
+
RayoDirecto = zeros(size(ListaArchivos,2),1);
35
+
RetFuente = zeros(size(ListaArchivos,2),1);
36
+
NumRayos = zeros(size(ListaArchivos,2),1);
37
+
Reverberacion = zeros(size(ListaArchivos,2),21);
38
+
Absorcion = zeros(size(ListaArchivos,2),21);
39
+
Directividad = zeros(size(ListaArchivos,2),21);
40
+
RuidoFondo = zeros(size(ListaArchivos,2),21);
41
+
SPLFuente = zeros(size(ListaArchivos,2),21);
42
+
handles.SPLm = cell(size(ListaArchivos,2),21);
43
+
for i=1:size(ListaArchivos,2)
44
+
45
+
% Distancia entre receptores y fuente
46
+
A = EASEPost2Distan(ListaArchivos{i});
47
+
Distancia(i) = sqrt(...
48
+
(A(1,1)-A(2,1))^2 + (A(1,2)-A(2,2))^2 + (A(1,3)-A(2,3))^2);
49
+
PosReceptor(i,:) = A(1,:);
50
+
PosFuente(i,:) = A(2,:);
51
+
52
+
% Nombre de la fuente y de los receptores
53
+
A = EASEPost2Info(ListaArchivos{i});
54
+
Fuente{i} = char(A(2));
55
+
Receptor{i} = char(A(1));
56
+
57
+
% Parametros de la sala
58
+
% ToDo: Almacenar la primera vez y el resto comprobar que son iguales, si
59
+
% no son iguales mostrar error y detener la importaci�n
60
+
A = EASEPost2Sala(ListaArchivos{i});
61
+
Volumen = A(1);
62
+
Superficie = A(2);
63
+
Temperatura = A(3);
64
+
Humedad = A(4);
65
+
66
+
% Datos varios
67
+
A = EASEPost2Otros(ListaArchivos{i});
68
+
RayoDirecto(i) = A(1);
69
+
RetFuente(i) = A(2);
70
+
NumRayos(i) = A(3);
71
+
72
+
% RT, Absorcion de la sala, Directividad fuente
73
+
A = EASEPost2Param(ListaArchivos{i});
74
+
Reverberacion(i,:) = A(1,:);
75
+
Absorcion(i,:) = A(2,:);
76
+
Directividad(i,:) = A(3,:);
77
+
78
+
% Ruido de fondo, SPL fuente
79
+
A = EASEPost2Param2(ListaArchivos{i});
80
+
RuidoFondo(i,:) = A(1,:);
81
+
SPLFuente(i,:) = A(2,:);
82
+
83
+
% Obtener historia temporal para cada receptor
84
+
% Obtiene una matriz con la columna de delay y los niveles en tercios de
85
+
% octava
86
+
A = EASEPost2Matrix(ListaArchivos{i});
87
+
% Eliminar el delay de llegada, restando el valor minimo de llegada a todo
88
+
% el vector de tiempos. Asi aseguramos que la integracion por rangos
89
+
% temporales sea correcta.
90
+
Tnew = A(:,1) - RayoDirecto(i);
91
+
% Almacenar en el cell 'SPLm' la informacion del archivo actual 'i'
92
+
handles.SPLm{i,1} = []; % Aqui ira el SPL por milisegundo
93
+
handles.SPLm{i,2} = Fuente{i};
94
+
handles.SPLm{i,3} = Receptor{i};
95
+
handles.SPLm{i,4} = PosFuente(i,:);
96
+
handles.SPLm{i,5} = PosReceptor(i,:);
97
+
handles.SPLm{i,6} = Distancia(i);
98
+
handles.SPLm{i,7} = RetFuente(i);
99
+
handles.SPLm{i,8} = NumRayos(i);
100
+
handles.SPLm{i,9} = RayoDirecto(i);
101
+
handles.SPLm{i,10} = Reverberacion(i,:);
102
+
handles.SPLm{i,11} = Absorcion(i,:);
103
+
handles.SPLm{i,12} = RuidoFondo(i,:);
104
+
handles.SPLm{i,13} = SPLFuente(i,:);
105
+
handles.SPLm{i,14} = Directividad(i,:);
106
+
handles.SPLm{i,15} = Volumen;
107
+
handles.SPLm{i,16} = Superficie;
108
+
handles.SPLm{i,17} = Temperatura;
109
+
handles.SPLm{i,18} = Humedad;
110
+
% Historia temporal completa
111
+
handles.SPLm{i,19} = [Tnew,A(:,2:end)];
112
+
% Historia temporal procesada por milisegundo en bandas de octava
113
+
handles.SPLm{i,20} = [];
114
+
% Historia temporal procesada por milisegundo en tercios de octava
115
+
handles.SPLm{i,21} = [];
116
+
117
+
end
118
+
119
+
% En primer lugar agrupa los niveles en bloques de 1 milisegundo, ya que hay
120
+
% varios niveles por milisegundo, y promedia los valores de cada bloque
121
+
for i=1:numel(handles.SPLm(:,1))
122
+
123
+
124
+
for j=1:max(ceil(handles.SPLm{i,19}(:,1)))
125
+
Ini = find(handles.SPLm{i,19}(:,1)>=j-1,1); % Inicio del bloque de 'j' milisegundos
126
+
Fin = find(handles.SPLm{i,19}(:,1)>j,1)-1; % Fin del bloque de 'j' milisegundos
127
+
% Si ha llegado al final del vector, fin estara vacio y inicio y fin
128
+
% seran iguales para obtener el valor ultimo del vector
129
+
if isempty(find(handles.SPLm{i,19}(:,1)>j,1)); Fin=Ini; end
130
+
% Si el bloque del milisegundo 'j' existe, ini sera menor o igual a fin,
131
+
% si no existe ini sera mayor a fin y no habra valor en ese milisegundo.
132
+
if Ini<=Fin
133
+
Valores = handles.SPLm{i,19}(Ini:Fin,2:end);
134
+
135
+
if max(Valores(:)) > 0
136
+
% Suma el bloque por octavas y vuelve a sumar para obtener un valor
137
+
ValMS(j) = 10*log10(sum(sum(10.^(Valores/10),1)));
138
+
% Suma los rayos recibidos en el milisegundo 'j', y lo pasa de
139
+
% tercios de octava a bandas de octava
140
+
ValMSthirdoct(j,:) = sum(10.^(Valores/10),1);
141
+
ValMSoct(j,:) = 10*log10([sum(ValMSthirdoct(j,1:3)),... % 125 Hz
142
+
sum(ValMSthirdoct(j,4:6)),... % 250 Hz
143
+
sum(ValMSthirdoct(j,7:9)),... % 500 Hz
144
+
sum(ValMSthirdoct(j,10:12)),...% 1 kHz
145
+
sum(ValMSthirdoct(j,13:15)),...% 2 kHz
146
+
sum(ValMSthirdoct(j,16:18)),...% 4 kHz
147
+
sum(ValMSthirdoct(j,19:21))]); % 8 kHz
148
+
ValMSthirdoct(j,:) = 10*log10(ValMSthirdoct(j,:));
149
+
else
150
+
ValMS(j) = 0;
151
+
ValMSoct(j,:) = 0;
152
+
ValMSthirdoct(j,:) = 0;
153
+
end
154
+
else
155
+
ValMS(j)=NaN;
156
+
ValMSoct(j,:) = NaN;
157
+
ValMSthirdoct(j,:) = NaN;
158
+
end
159
+
160
+
end
161
+
% Todos los valores 0 y menores a 0 se dejan vacios
162
+
%ValMS(ValMS==0) = NaN;
163
+
%ValMSoct(ValMSoct==0) = NaN;
164
+
% Si se utiliza matlab 2014b o anterior realiza la interpolacion
165
+
% manualmente, sino utiliza la funcion fillmising
166
+
if verLessThan('matlab','8.7')
167
+
% Interpolacion manual global
168
+
bd=isnan(ValMS);
169
+
gd=find(~bd);
170
+
bd([1:(min(gd)-1) (max(gd)+1):end])=0;
171
+
ValMS(bd)=interp1(gd,ValMS(gd),find(bd));
172
+
% Interpolacion manual por bandas
173
+
for jj=1:7
174
+
bd=isnan(ValMSoct(:,jj));
175
+
gd=find(~bd);
176
+
bd([1:(min(gd)-1) (max(gd)+1):end])=0;
177
+
ValMSoct(bd,jj)=interp1(gd,ValMSoct(gd,jj),find(bd));
178
+
end
179
+
for jj=1:21
180
+
bd=isnan(ValMSthirdoct(:,jj));
181
+
gd=find(~bd);
182
+
bd([1:(min(gd)-1) (max(gd)+1):end])=0;
183
+
ValMSthirdoct(bd,jj)=interp1(gd,ValMSthirdoct(gd,jj),find(bd));
184
+
end
185
+
else
186
+
% Interpolar valores en los milisegundos que haya valor 0, desde el inicio
187
+
% al final del vector
188
+
ValMS = fillmissing(ValMS,'linear','EndValues','extrap');
189
+
ValMSoct = fillmissing(ValMSoct,'linear','EndValues','extrap');
190
+
ValMSthirdoct = fillmissing(ValMSthirdoct,'linear','EndValues','extrap');
191
+
end
192
+
handles.SPLm{i,1} = ValMS';
193
+
handles.SPLm{i,20} = ValMSoct;
194
+
handles.SPLm{i,21} = ValMSthirdoct;
195
+
end
196
+
% Obtener el tiempo maximo valor elegible para representacion. Obtiene el
197
+
% maximo de cada receptor y despues busca el minimo entre estos.
198
+
handles.MaxElegible = min(cell2mat(cellfun( @(val) numel(val), handles.SPLm(:,1),'UniformOutput',0)))-1;
199
+
200
+
% Ordenar los datos segun la ID de fuente y despues por receptor. Asi se
201
+
% tiene agrupados los datos segun la fuente, necesario para elegir una
202
+
% fuente u otra en el programa.
203
+
handles.SPLm = sortrows(handles.SPLm,[2,3]);
204
+
% Actualizar variable en la aplicacion
205
+
guidata(hObject, handles);
206
+
207
+
end
208
+
+824
include/multiWaitbar.m
+824
include/multiWaitbar.m
···
1
+
function cancel = multiWaitbar( handles,label, varargin )
2
+
%multiWaitbar: add, remove or update an entry on the multi waitbar
3
+
%
4
+
% multiWaitbar(LABEL,VALUE) adds a waitbar for the specified label, or
5
+
% if it already exists updates the value. LABEL must be a string and
6
+
% VALUE a number between zero and one or the string 'Close' to remove the
7
+
% entry Setting value equal to 0 or 'Reset' will cause the progress bar
8
+
% to reset and the time estimate to be re-initialized.
9
+
%
10
+
% multiWaitbar(LABEL,COMMAND,VALUE,...) or
11
+
% multiWaitbar(LABEL,VALUE,COMMAND,VALUE,...)
12
+
% passes one or more command/value pairs for changing the named waitbar
13
+
% entry. Possible commands include:
14
+
% 'Value' Set the value of the named waitbar entry. The
15
+
% corresponding value must be a number between 0 and 1.
16
+
% 'Increment' Increment the value of the named waitbar entry. The
17
+
% corresponding value must be a number between 0 and 1.
18
+
% 'Color' Change the color of the named waitbar entry. The
19
+
% value must be an RGB triple, e.g. [0.1 0.2 0.3], or a
20
+
% single-character color name, e.g. 'r', 'b', 'm'.
21
+
% 'Relabel' Change the label of the named waitbar entry. The
22
+
% value must be the new name.
23
+
% 'Reset' Set the named waitbar entry back to zero and reset its
24
+
% timer. No value need be specified.
25
+
% 'CanCancel' [on|off] should a "cancel" button be shown for this bar
26
+
% (default 'off').
27
+
% 'CancelFcn' Function to call in the event that the user cancels.
28
+
% 'ResetCancel' Reset the "cancelled" flag for an entry (ie. if you
29
+
% decide not to cancel).
30
+
% 'Close' Remove the named waitbar entry.
31
+
% 'Busy' Puts this waitbar in "busy mode" where a small bar
32
+
% bounces back and forth. Return to normal progress display
33
+
% using the 'Reset' command.
34
+
%
35
+
% cancel = multiWaitbar(LABEL,VALUE) also returns whether the user has
36
+
% clicked the "cancel" button for this entry (true or false). Two
37
+
% mechanisms are provided for cancelling an entry if the 'CanCancel'
38
+
% setting is 'on'. The first is just to check the return argument and if
39
+
% it is true abort the task. The second is to set a 'CancelFcn' that is
40
+
% called when the user clicks the cancel button, much as is done for
41
+
% MATLAB's built-in WAITBAR. In either case, you can use the
42
+
% 'ResetCancel' command if you don't want to cancel after all.
43
+
%
44
+
% multiWaitbar('CLOSEALL') closes the waitbar window.
45
+
%
46
+
% Example:
47
+
% multiWaitbar( 'CloseAll' );
48
+
% multiWaitbar( 'Task 1', 0 );
49
+
% multiWaitbar( 'Task 2', 0.5, 'Color', 'b' );
50
+
% multiWaitbar( 'Task 3', 'Busy');
51
+
% multiWaitbar( 'Task 1', 'Value', 0.1 );
52
+
% multiWaitbar( 'Task 2', 'Increment', 0.2 );
53
+
% multiWaitbar( 'Task 3', 'Reset' ); % Disables "busy" mode
54
+
% multiWaitbar( 'Task 3', 'Value', 0.3 );
55
+
% multiWaitbar( 'Task 2', 'Close' );
56
+
% multiWaitbar( 'Task 3', 'Close' );
57
+
% multiWaitbar( 'Task 1', 'Close' );
58
+
%
59
+
% Example:
60
+
% multiWaitbar( 'Task 1', 0, 'CancelFcn', @(a,b) disp( ['Cancel ',a] ) );
61
+
% for ii=1:100
62
+
% abort = multiWaitbar( 'Task 1', ii/100 );
63
+
% if abort
64
+
% % Here we would normally ask the user if they're sure
65
+
% break
66
+
% else
67
+
% pause( 1 )
68
+
% end
69
+
% end
70
+
% multiWaitbar( 'Task 1', 'Close' )
71
+
%
72
+
% Example:
73
+
% multiWaitbar( 'CloseAll' );
74
+
% multiWaitbar( 'Red...', 7/7, 'Color', [0.8 0.0 0.1] );
75
+
% multiWaitbar( 'Orange...', 6/7, 'Color', [1.0 0.4 0.0] );
76
+
% multiWaitbar( 'Yellow...', 5/7, 'Color', [0.9 0.8 0.2] );
77
+
% multiWaitbar( 'Green...', 4/7, 'Color', [0.2 0.9 0.3] );
78
+
% multiWaitbar( 'Blue...', 3/7, 'Color', [0.1 0.5 0.8] );
79
+
% multiWaitbar( 'Indigo...', 2/7, 'Color', [0.4 0.1 0.5] );
80
+
% multiWaitbar( 'Violet...', 1/7, 'Color', [0.8 0.4 0.9] );
81
+
%
82
+
% Thanks to Jesse Hopkins for suggesting the "busy" mode.
83
+
84
+
% Author: Ben Tordoff
85
+
% Copyright 2007-2014 The MathWorks, Inc.
86
+
87
+
% Modificaciones para CATT2Matlab
88
+
global titleprogress
89
+
global timedays
90
+
global timehours
91
+
global timemins
92
+
global timesecs
93
+
global timesec
94
+
titleprogress = handles.LPROGRESSTITLE;
95
+
timedays = handles.LTIMEDAYS;
96
+
timehours = handles.LTIMEHOURS;
97
+
timemins = handles.LTIMEMINS;
98
+
timesecs = handles.LTIMESECS;
99
+
timesec = handles.LTIMESEC;
100
+
101
+
102
+
persistent FIGH;
103
+
cancel = false;
104
+
105
+
% Check basic inputs
106
+
narginchk( 1, inf); % - kept for backwards compatibility
107
+
if ~ischar( label )
108
+
error( 'multiWaitbar:BadArg', 'LABEL must be the name of the progress entry (i.e. a string)' );
109
+
end
110
+
111
+
% Try to get hold of the figure
112
+
if isempty( FIGH ) || ~ishandle( FIGH )
113
+
FIGH = findall( 0, 'Type', 'figure', 'Tag', 'multiWaitbar:Figure' );
114
+
if isempty(FIGH)
115
+
FIGH = iCreateFig();
116
+
else
117
+
FIGH = handle( FIGH(1) );
118
+
end
119
+
end
120
+
121
+
% Check for close all and stop early
122
+
if any( strcmpi( label, {'CLOSEALL','CLOSE ALL'} ) )
123
+
iDeleteFigure(FIGH);
124
+
return;
125
+
end
126
+
127
+
% Make sure we're on-screen
128
+
if ~strcmpi( FIGH.Visible, 'on' )
129
+
FIGH.Visible = 'on';
130
+
end
131
+
132
+
% Get the list of entries and see if this one already exists
133
+
entries = getappdata( FIGH, 'ProgressEntries' );
134
+
if isempty(entries)
135
+
idx = [];
136
+
else
137
+
idx = find( strcmp( label, {entries.Label} ), 1, 'first' );
138
+
end
139
+
bgcol = getappdata( FIGH, 'DefaultProgressBarBackgroundColor' );
140
+
141
+
% If it doesn't exist, create it
142
+
needs_redraw = false;
143
+
entry_added = isempty(idx);
144
+
if entry_added
145
+
% Create a new entry
146
+
defbarcolor = getappdata( FIGH, 'DefaultProgressBarColor' );
147
+
entries = iAddEntry( FIGH, entries, label, 0, defbarcolor, bgcol );
148
+
idx = numel( entries );
149
+
end
150
+
151
+
% Check if the user requested a cancel
152
+
if nargout
153
+
cancel = entries(idx).Cancel;
154
+
end
155
+
156
+
% Parse the inputs. We shortcut the most common case as an efficiency
157
+
force_update = false;
158
+
if nargin==2 && isnumeric( varargin{1} )
159
+
entries(idx).LastValue = entries(idx).Value;
160
+
entries(idx).Value = max( 0, min( 1, varargin{1} ) );
161
+
entries(idx).Busy = false;
162
+
needs_update = true;
163
+
else
164
+
[params,values] = iParseInputs( varargin{:} );
165
+
166
+
needs_update = false;
167
+
for ii=1:numel( params )
168
+
switch upper( params{ii} )
169
+
case 'BUSY'
170
+
entries(idx).Busy = true;
171
+
needs_update = true;
172
+
173
+
case 'VALUE'
174
+
entries(idx).LastValue = entries(idx).Value;
175
+
entries(idx).Value = max( 0, min( 1, values{ii} ) );
176
+
entries(idx).Busy = false;
177
+
needs_update = true;
178
+
179
+
case {'INC','INCREMENT'}
180
+
entries(idx).LastValue = entries(idx).Value;
181
+
entries(idx).Value = max( 0, min( 1, entries(idx).Value + values{ii} ) );
182
+
entries(idx).Busy = false;
183
+
needs_update = true;
184
+
185
+
case {'COLOR','COLOUR'}
186
+
entries(idx).CData = iMakeColors( values{ii}, 16 );
187
+
needs_update = true;
188
+
force_update = true;
189
+
190
+
case {'RELABEL', 'UPDATELABEL'}
191
+
% Make sure we have a string as the value and that it
192
+
% doesn't already appear
193
+
if ~ischar( values{ii} )
194
+
error( 'multiWaitbar:BadString', 'Value for ''Relabel'' must be a string.' );
195
+
end
196
+
if ismember( values{ii}, {entries.Label} )
197
+
error( 'multiWaitbar:NameAlreadyExists', 'Cannot relabel an entry to a label that already exists.' );
198
+
end
199
+
entries(idx).Label = values{ii};
200
+
needs_update = true;
201
+
force_update = true;
202
+
203
+
case {'CANCANCEL'}
204
+
if ~ischar( values{ii} ) || ~any( strcmpi( values{ii}, {'on','off'} ) )
205
+
error( 'multiWaitbar:BadString', 'Parameter ''CanCancel'' must be a ''on'' or ''off''.' );
206
+
end
207
+
entries(idx).CanCancel = strcmpi( values{ii}, 'on' );
208
+
entries(idx).Cancel = false;
209
+
needs_redraw = true;
210
+
211
+
case {'RESETCANCEL'}
212
+
entries(idx).Cancel = false;
213
+
needs_redraw = true;
214
+
215
+
case {'CANCELFCN'}
216
+
if ~isa( values{ii}, 'function_handle' )
217
+
error( 'multiWaitbar:BadFunction', 'Parameter ''CancelFcn'' must be a valid function handle.' );
218
+
end
219
+
entries(idx).CancelFcn = values{ii};
220
+
if ~entries(idx).CanCancel
221
+
entries(idx).CanCancel = true;
222
+
end
223
+
needs_redraw = true;
224
+
225
+
case {'CLOSE','DONE'}
226
+
if ~isempty(idx)
227
+
% Remove the selected entry
228
+
entries = iDeleteEntry( entries, idx );
229
+
end
230
+
if isempty( entries )
231
+
iDeleteFigure( FIGH );
232
+
% With the window closed, there's nothing else to do
233
+
return;
234
+
else
235
+
needs_redraw = true;
236
+
end
237
+
% We can't continue after clearing the entry, so jump out
238
+
break;
239
+
240
+
otherwise
241
+
error( 'multiWaitbar:BadArg', 'Unrecognized command: ''%s''', params{ii} );
242
+
243
+
end
244
+
end
245
+
end
246
+
247
+
% Now work out what to update/redraw
248
+
if needs_redraw
249
+
setappdata( FIGH, 'ProgressEntries', entries );
250
+
iRedraw( FIGH );
251
+
% NB: Redraw includes updating all bars, so never need to do both
252
+
elseif needs_update
253
+
[entries(idx),needs_redraw] = iUpdateEntry( entries(idx), force_update );
254
+
setappdata( FIGH, 'ProgressEntries', entries );
255
+
% NB: if anything was updated onscreen, "needs_redraw" is now true.
256
+
end
257
+
if entry_added || needs_redraw
258
+
% If the shape or size has changed, do a full redraw, including events
259
+
drawnow();
260
+
end
261
+
262
+
% If we have any "busy" entries, start the timer, otherwise stop it.
263
+
myTimer = getappdata( FIGH, 'BusyTimer' );
264
+
if any([entries.Busy])
265
+
if strcmpi(myTimer.Running,'off')
266
+
start(myTimer);
267
+
end
268
+
else
269
+
if strcmpi(myTimer.Running,'on')
270
+
stop(myTimer);
271
+
end
272
+
end
273
+
274
+
clearvars -global titleprogress timedays timehours timemins timesecs timesec
275
+
end % multiWaitbar
276
+
277
+
278
+
%-------------------------------------------------------------------------%
279
+
function [params, values] = iParseInputs( varargin )
280
+
% Parse the input arguments, extracting a list of commands and values
281
+
idx = 1;
282
+
params = {};
283
+
values = {};
284
+
if nargin==0
285
+
return;
286
+
end
287
+
if isnumeric( varargin{1} )
288
+
params{idx} = 'Value';
289
+
values{idx} = varargin{1};
290
+
idx = idx + 1;
291
+
end
292
+
293
+
while idx <= nargin
294
+
param = varargin{idx};
295
+
if ~ischar( param )
296
+
error( 'multiWaitbar:BadSyntax', 'Additional properties must be supplied as property-value pairs' );
297
+
end
298
+
params{end+1,1} = param; %#ok<AGROW>
299
+
values{end+1,1} = []; %#ok<AGROW>
300
+
switch upper( param )
301
+
case {'DONE','CLOSE','RESETCANCEL'}
302
+
% No value needed, and stop
303
+
break;
304
+
case {'BUSY'}
305
+
% No value needed but parsing should continue
306
+
idx = idx + 1;
307
+
case {'RESET','ZERO','SHOW'}
308
+
% All equivalent to saying ('Value', 0)
309
+
params{end} = 'Value';
310
+
values{end} = 0;
311
+
idx = idx + 1;
312
+
otherwise
313
+
if idx==nargin
314
+
error( 'multiWaitbar:BadSyntax', 'Additional properties must be supplied as property-value pairs' );
315
+
end
316
+
values{end,1} = varargin{idx+1};
317
+
idx = idx + 2;
318
+
end
319
+
end
320
+
if isempty( params )
321
+
error( 'multiWaitbar:BadSyntax', 'Must specify a value or a command' );
322
+
end
323
+
end % iParseInputs
324
+
325
+
%-------------------------------------------------------------------------%
326
+
function fobj = iCreateFig()
327
+
global titleprogress
328
+
% Create the progress bar group window
329
+
bgcol = get(0,'DefaultUIControlBackgroundColor');
330
+
f = figure( ...
331
+
'Name', titleprogress, ...
332
+
'Tag', 'multiWaitbar:Figure', ...
333
+
'Color', bgcol, ...
334
+
'MenuBar', 'none', ...
335
+
'ToolBar', 'none', ...
336
+
'WindowStyle', 'normal', ... % We don't want to be docked!
337
+
'HandleVisibility', 'off', ...
338
+
'IntegerHandle', 'off', ...
339
+
'Visible', 'off', ...
340
+
'NumberTitle', 'off' );
341
+
% Resize and centre on the first screen
342
+
screenSize = get(0,'ScreenSize');
343
+
figSz = [360 42];
344
+
figPos = ceil((screenSize(1,3:4)-figSz)/2);
345
+
fobj = handle( f );
346
+
fobj.Position = [figPos, figSz];
347
+
setappdata( fobj, 'ProgressEntries', [] );
348
+
% Make sure we have the image
349
+
defbarcolor = [0.8 0.0 0.1];
350
+
barbgcol = uint8( 255*0.75*bgcol );
351
+
setappdata( fobj, 'DefaultProgressBarBackgroundColor', barbgcol );
352
+
setappdata( fobj, 'DefaultProgressBarColor', defbarcolor );
353
+
setappdata( fobj, 'DefaultProgressBarSize', [350 16] );
354
+
% Create the timer to use for "Busy" mode, being sure to delete any
355
+
% existing ones
356
+
delete( timerfind('Tag', 'MultiWaitbarTimer') );
357
+
myTimer = timer( ...
358
+
'TimerFcn', @(src,evt) iTimerFcn(f), ...
359
+
'Period', 0.02, ...
360
+
'ExecutionMode', 'FixedRate', ...
361
+
'Tag', 'MultiWaitbarTimer' );
362
+
setappdata( fobj, 'BusyTimer', myTimer );
363
+
364
+
% Setup the resize function after we've finished setting up the figure to
365
+
% avoid excessive redraws
366
+
fobj.ResizeFcn = @iRedraw;
367
+
fobj.CloseRequestFcn = @iCloseFigure;
368
+
end % iCreateFig
369
+
370
+
%-------------------------------------------------------------------------%
371
+
function cdata = iMakeColors( baseColor, height )
372
+
% Creates a shiny bar from a single base color
373
+
lightColor = [1 1 1];
374
+
badColorErrorID = 'multiWaitbar:BadColor';
375
+
badColorErrorMsg = 'Colors must be a three element vector [R G B] or a single character (''r'', ''g'' etc.)';
376
+
377
+
if ischar(baseColor)
378
+
switch upper(baseColor)
379
+
case 'K'
380
+
baseColor = [0.1 0.1 0.1];
381
+
case 'R'
382
+
baseColor = [0.8 0 0];
383
+
case 'G'
384
+
baseColor = [0 0.6 0];
385
+
case 'B'
386
+
baseColor = [0 0 0.8];
387
+
case 'C'
388
+
baseColor = [0.2 0.8 0.9];
389
+
case 'M'
390
+
baseColor = [0.6 0 0.6];
391
+
case 'Y'
392
+
baseColor = [0.9 0.8 0.2];
393
+
case 'W'
394
+
baseColor = [0.9 0.9 0.9];
395
+
otherwise
396
+
error( badColorErrorID, badColorErrorMsg );
397
+
end
398
+
else
399
+
if numel(baseColor) ~= 3
400
+
error( badColorErrorID, badColorErrorMsg );
401
+
end
402
+
if isa( baseColor, 'uint8' )
403
+
baseColor = double( baseColor ) / 255;
404
+
elseif isa( baseColor, 'double' )
405
+
if any(baseColor>1) || any(baseColor<0)
406
+
error( 'multiWaitbar:BadColorValue', 'Color values must be in the range 0 to 1 inclusive.' );
407
+
end
408
+
else
409
+
error( badColorErrorID, badColorErrorMsg );
410
+
end
411
+
end
412
+
413
+
% By this point we should have a double precision 3-element vector.
414
+
cols = repmat( baseColor, [height, 1] );
415
+
416
+
breaks = max( 1, round( height * [1 25 50 75 88 100] / 100 ) );
417
+
cols(breaks(1),:) = 0.6*baseColor;
418
+
cols(breaks(2),:) = lightColor - 0.4*(lightColor-baseColor);
419
+
cols(breaks(3),:) = baseColor;
420
+
cols(breaks(4),:) = min( baseColor*1.2, 1.0 );
421
+
cols(breaks(5),:) = min( baseColor*1.4, 0.95 ) + 0.05;
422
+
cols(breaks(6),:) = min( baseColor*1.6, 0.9 ) + 0.1;
423
+
424
+
y = 1:height;
425
+
cols(:,1) = max( 0, min( 1, interp1( breaks, cols(breaks,1), y, 'pchip' ) ) );
426
+
cols(:,2) = max( 0, min( 1, interp1( breaks, cols(breaks,2), y, 'pchip' ) ) );
427
+
cols(:,3) = max( 0, min( 1, interp1( breaks, cols(breaks,3), y, 'pchip' ) ) );
428
+
cdata = uint8( 255 * cat( 3, cols(:,1), cols(:,2), cols(:,3) ) );
429
+
end % iMakeColors
430
+
431
+
432
+
%-------------------------------------------------------------------------%
433
+
function cdata = iMakeBackground( baseColor, height )
434
+
% Creates a shaded background
435
+
if isa( baseColor, 'uint8' )
436
+
baseColor = double( baseColor ) / 255;
437
+
end
438
+
439
+
ratio = 1 - exp( -0.5-2*(1:height)/height )';
440
+
cdata = uint8( 255 * cat( 3, baseColor(1)*ratio, baseColor(2)*ratio, baseColor(3)*ratio ) );
441
+
end % iMakeBackground
442
+
443
+
%-------------------------------------------------------------------------%
444
+
function entries = iAddEntry( parent, entries, label, value, color, bgcolor )
445
+
% Add a new entry to the progress bar
446
+
447
+
% Create bar coloring
448
+
psize = getappdata( parent, 'DefaultProgressBarSize' );
449
+
cdata = iMakeColors( color, 16 );
450
+
% Create background image
451
+
barcdata = iMakeBackground( bgcolor, psize(2) );
452
+
453
+
% Work out the size in advance
454
+
labeltext = uicontrol( 'Style', 'Text', ...
455
+
'String', label, ...
456
+
'Parent', parent, ...
457
+
'HorizontalAlignment', 'Left' );
458
+
etatext = uicontrol( 'Style', 'Text', ...
459
+
'String', '', ...
460
+
'Parent', parent, ...
461
+
'HorizontalAlignment', 'Right' );
462
+
progresswidget = uicontrol( 'Style', 'Checkbox', ...
463
+
'String', '', ...
464
+
'Parent', parent, ...
465
+
'Position', [5 5 psize], ...
466
+
'CData', barcdata );
467
+
cancelwidget = uicontrol( 'Style', 'PushButton', ...
468
+
'String', '', ...
469
+
'FontWeight', 'Bold', ...
470
+
'Parent', parent, ...
471
+
'Position', [5 5 16 16], ...
472
+
'CData', iMakeCross( 8 ), ...
473
+
'Callback', @(src,evt) iCancelEntry( src, label ), ...
474
+
'Visible', 'off' );
475
+
mypanel = uipanel( 'Parent', parent, 'Units', 'Pixels' );
476
+
477
+
newentry = struct( ...
478
+
'Label', label, ...
479
+
'Value', value, ...
480
+
'LastValue', inf, ...
481
+
'Created', tic(), ...
482
+
'LabelText', labeltext, ...
483
+
'ETAText', etatext, ...
484
+
'ETAString', '', ...
485
+
'Progress', progresswidget, ...
486
+
'ProgressSize', psize, ...
487
+
'Panel', mypanel, ...
488
+
'BarCData', barcdata, ...
489
+
'CData', cdata, ...
490
+
'BackgroundCData', barcdata, ...
491
+
'CanCancel', false, ...
492
+
'CancelFcn', [], ...
493
+
'CancelButton', cancelwidget, ...
494
+
'Cancel', false, ...
495
+
'Busy', false );
496
+
if isempty( entries )
497
+
entries = newentry;
498
+
else
499
+
entries = [entries;newentry];
500
+
end
501
+
% Store in figure before the redraw
502
+
setappdata( parent, 'ProgressEntries', entries );
503
+
if strcmpi( get( parent, 'Visible' ), 'on' )
504
+
iRedraw( parent, [] );
505
+
else
506
+
set( parent, 'Visible', 'on' );
507
+
end
508
+
end % iAddEntry
509
+
510
+
%-------------------------------------------------------------------------%
511
+
function entries = iDeleteEntry( entries, idx )
512
+
delete( entries(idx).LabelText );
513
+
delete( entries(idx).ETAText );
514
+
delete( entries(idx).CancelButton );
515
+
delete( entries(idx).Progress );
516
+
delete( entries(idx).Panel );
517
+
entries(idx,:) = [];
518
+
end % iDeleteEntry
519
+
520
+
%-------------------------------------------------------------------------%
521
+
function entries = iCancelEntry( src, name )
522
+
figh = ancestor( src, 'figure' );
523
+
entries = getappdata( figh, 'ProgressEntries' );
524
+
if isempty(entries)
525
+
% The entries have been lost - nothing can be done.
526
+
return
527
+
end
528
+
idx = find( strcmp( name, {entries.Label} ), 1, 'first' );
529
+
530
+
% Set the cancel flag so that the user is told on next update
531
+
entries(idx).Cancel = true;
532
+
setappdata( figh, 'ProgressEntries', entries );
533
+
534
+
% If a user function is supplied, call it
535
+
if ~isempty( entries(idx).CancelFcn )
536
+
feval( entries(idx).CancelFcn, name, 'Cancelled' );
537
+
end
538
+
539
+
end % iCancelEntry
540
+
541
+
542
+
%-------------------------------------------------------------------------%
543
+
function [entry,updated] = iUpdateEntry( entry, force )
544
+
% Update one progress bar
545
+
546
+
% Deal with busy entries separately
547
+
if entry.Busy
548
+
entry = iUpdateBusyEntry(entry);
549
+
updated = true;
550
+
return;
551
+
end
552
+
553
+
% Some constants
554
+
marker_weight = 0.8;
555
+
556
+
% Check if the label needs updating
557
+
updated = force;
558
+
val = entry.Value;
559
+
lastval = entry.LastValue;
560
+
561
+
% Now update the bar
562
+
psize = entry.ProgressSize;
563
+
filled = max( 1, round( val*psize(1) ) );
564
+
lastfilled = max( 1, round( lastval*psize(1) ) );
565
+
566
+
% We do some careful checking so that we only redraw what we have to. This
567
+
% makes a small speed difference, but every little helps!
568
+
if force || (filled<lastfilled)
569
+
% Create the bar background
570
+
startIdx = 1;
571
+
bgim = entry.BackgroundCData(:,ones( 1, psize(1)-filled ),:);
572
+
barim = iMakeBarImage(entry.CData, startIdx, filled);
573
+
progresscdata = [barim,bgim];
574
+
575
+
% Add light/shadow around the markers
576
+
markers = round( (0.1:0.1:val)*psize(1) );
577
+
markers(markers<startIdx | markers>(filled-2)) = [];
578
+
highlight = [marker_weight*entry.CData, 255 - marker_weight*(255-entry.CData)];
579
+
for ii=1:numel( markers )
580
+
progresscdata(:,markers(ii)+[-1,0],:) = highlight;
581
+
end
582
+
583
+
% Set the image into the checkbox
584
+
entry.BarCData = progresscdata;
585
+
set( entry.Progress, 'cdata', progresscdata );
586
+
updated = true;
587
+
588
+
elseif filled > lastfilled
589
+
% Just need to update the existing data
590
+
progresscdata = entry.BarCData;
591
+
startIdx = max(1,lastfilled-1);
592
+
% Repmat is the obvious way to fill the bar, but BSXFUN is often
593
+
% faster. Indexing is obscure but faster still.
594
+
progresscdata(:,startIdx:filled,:) = iMakeBarImage(entry.CData, startIdx, filled);
595
+
596
+
% Add light/shadow around the markers
597
+
markers = round( (0.1:0.1:val)*psize(1) );
598
+
markers(markers<startIdx | markers>(filled-2)) = [];
599
+
highlight = [marker_weight*entry.CData, 255 - marker_weight*(255-entry.CData)];
600
+
for ii=1:numel( markers )
601
+
progresscdata(:,markers(ii)+[-1,0],:) = highlight;
602
+
end
603
+
604
+
entry.BarCData = progresscdata;
605
+
set( entry.Progress, 'CData', progresscdata );
606
+
updated = true;
607
+
end
608
+
609
+
% As an optimization, don't update any text if the bar didn't move and the
610
+
% percentage hasn't changed
611
+
decval = round( val*100 );
612
+
lastdecval = round( lastval*100 );
613
+
614
+
if ~updated && (decval == lastdecval)
615
+
return
616
+
end
617
+
618
+
% Now work out the remaining time
619
+
minTime = 3; % secs
620
+
if val <= 0
621
+
% Zero value, so clear the eta
622
+
entry.Created = tic();
623
+
elapsedtime = 0;
624
+
etaString = '';
625
+
else
626
+
elapsedtime = round(toc( entry.Created )); % in seconds
627
+
628
+
% Only show the remaining time if we've had time to estimate
629
+
if elapsedtime < minTime
630
+
% Not enough time has passed since starting, so leave blank
631
+
etaString = '';
632
+
else
633
+
% Calculate a rough ETA
634
+
eta = elapsedtime * (1-val) / val;
635
+
etaString = iGetTimeString( eta );
636
+
end
637
+
end
638
+
639
+
if ~isequal( etaString, entry.ETAString )
640
+
set( entry.ETAText, 'String', etaString );
641
+
entry.ETAString = etaString;
642
+
updated = true;
643
+
end
644
+
645
+
% Update the label too
646
+
if force || elapsedtime > minTime
647
+
if force || (decval ~= lastdecval)
648
+
labelstr = [entry.Label, sprintf( ' (%d%%)', decval )];
649
+
set( entry.LabelText, 'String', labelstr );
650
+
updated = true;
651
+
end
652
+
end
653
+
654
+
end % iUpdateEntry
655
+
656
+
function eta = iGetTimeString( remainingtime )
657
+
global timedays
658
+
global timehours
659
+
global timemins
660
+
global timesecs
661
+
global timesec
662
+
if remainingtime > 172800 % 2 days
663
+
eta = [sprintf( '%d ', round(remainingtime/86400) ),timedays];
664
+
else
665
+
if remainingtime > 7200 % 2 hours
666
+
eta = [sprintf( '%d ', round(remainingtime/3600) ),timehours];
667
+
else
668
+
if remainingtime > 120 % 2 mins
669
+
eta = [sprintf( '%d ', round(remainingtime/60) ),timemins];
670
+
else
671
+
% Seconds
672
+
remainingtime = round( remainingtime );
673
+
if remainingtime > 1
674
+
eta = [sprintf( '%d ', remainingtime ),timesecs];
675
+
elseif remainingtime == 1
676
+
eta = ['1 ',timesec];
677
+
else
678
+
eta = ''; % Nearly done (<1sec)
679
+
end
680
+
end
681
+
end
682
+
end
683
+
end % iGetTimeString
684
+
685
+
686
+
%-------------------------------------------------------------------------%
687
+
function entry = iUpdateBusyEntry( entry )
688
+
% Update a "busy" progress bar
689
+
% Make sure the widget is still OK
690
+
if ~ishandle(entry.Progress)
691
+
return
692
+
end
693
+
% Work out the new position. Since the bar is 0.1 long and needs to bounce,
694
+
% the position varies from 0 up to 0.9 then back down again. We achieve
695
+
% this with judicious use of "mod" with 1.8.
696
+
entry.Value = mod(entry.Value+0.01,1.8);
697
+
val = entry.Value;
698
+
if val>0.9
699
+
% Moving backwards
700
+
val = 1.8-val;
701
+
end
702
+
psize = entry.ProgressSize;
703
+
startIdx = max( 1, round( val*psize(1) ) );
704
+
endIdx = max( 1, round( (val+0.1)*psize(1) ) );
705
+
barLength = endIdx - startIdx + 1;
706
+
707
+
% Create the image
708
+
bgim = entry.BackgroundCData(:,ones( 1, psize(1) ),:);
709
+
barim = iMakeBarImage(entry.CData, 1, barLength);
710
+
bgim(:,startIdx:endIdx,:) = barim;
711
+
712
+
% Put it into the widget
713
+
entry.BarCData = bgim;
714
+
set( entry.Progress, 'CData', bgim );
715
+
end % iUpdateBusyEntry
716
+
717
+
718
+
%-------------------------------------------------------------------------%
719
+
function barim = iMakeBarImage(strip, startIdx, endIdx)
720
+
shadow1_weight = 0.4;
721
+
shadow2_weight = 0.7;
722
+
barLength = endIdx - startIdx + 1;
723
+
% Repmat is the obvious way to fill the bar, but BSXFUN is often
724
+
% faster. Indexing is obscure but faster still.
725
+
barim = strip(:,ones(1, barLength),:);
726
+
% Add highlight to the start of the bar
727
+
if startIdx <= 2 && barLength>=2
728
+
barim(:,1,:) = 255 - shadow1_weight*(255-strip);
729
+
barim(:,2,:) = 255 - shadow2_weight*(255-strip);
730
+
end
731
+
% Add shadow to the end of the bar
732
+
if endIdx>=4 && barLength>=2
733
+
barim(:,end,:) = shadow1_weight*strip;
734
+
barim(:,end-1,:) = shadow2_weight*strip;
735
+
end
736
+
end % iMakeBarImage
737
+
738
+
%-------------------------------------------------------------------------%
739
+
function iCloseFigure( fig, evt ) %#ok<INUSD>
740
+
% Closing the figure just makes it invisible
741
+
set( fig, 'Visible', 'off' );
742
+
end % iCloseFigure
743
+
744
+
%-------------------------------------------------------------------------%
745
+
function iDeleteFigure( fig )
746
+
% Actually destroy the figure
747
+
busyTimer = getappdata( fig, 'BusyTimer' );
748
+
stop( busyTimer );
749
+
delete( busyTimer );
750
+
delete( fig );
751
+
end % iDeleteFigure
752
+
753
+
%-------------------------------------------------------------------------%
754
+
function iRedraw( fig, evt ) %#ok<INUSD>
755
+
entries = getappdata( fig, 'ProgressEntries' );
756
+
fobj = handle( fig );
757
+
p = fobj.Position;
758
+
% p = get( fig, 'Position' );
759
+
border = 5;
760
+
textheight = 16;
761
+
barheight = 16;
762
+
panelheight = 10;
763
+
N = max( 1, numel( entries ) );
764
+
765
+
% Check the height is correct
766
+
heightperentry = textheight+barheight+panelheight;
767
+
requiredheight = 2*border + N*heightperentry - panelheight;
768
+
if ~isequal( p(4), requiredheight )
769
+
p(2) = p(2) + p(4) - requiredheight;
770
+
p(4) = requiredheight;
771
+
% In theory setting the position should re-trigger this callback, but
772
+
% in practice it doesn't, probably because we aren't calling "drawnow".
773
+
set( fig, 'Position', p )
774
+
end
775
+
ypos = p(4) - border;
776
+
width = p(3) - 2*border;
777
+
setappdata( fig, 'DefaultProgressBarSize', [width barheight] );
778
+
779
+
for ii=1:numel( entries )
780
+
set( entries(ii).LabelText, 'Position', [border ypos-textheight width*0.75 textheight] );
781
+
set( entries(ii).ETAText, 'Position', [border+width*0.75 ypos-textheight width*0.25 textheight] );
782
+
ypos = ypos - textheight;
783
+
if entries(ii).CanCancel
784
+
set( entries(ii).Progress, 'Position', [border ypos-barheight width-barheight+1 barheight] );
785
+
entries(ii).ProgressSize = [width-barheight barheight];
786
+
set( entries(ii).CancelButton, 'Visible', 'on', 'Position', [p(3)-border-barheight ypos-barheight barheight barheight] );
787
+
else
788
+
set( entries(ii).Progress, 'Position', [border ypos-barheight width+1 barheight] );
789
+
entries(ii).ProgressSize = [width barheight];
790
+
set( entries(ii).CancelButton, 'Visible', 'off' );
791
+
end
792
+
ypos = ypos - barheight;
793
+
set( entries(ii).Panel, 'Position', [-500 ypos-500-panelheight/2 p(3)+1000 500] );
794
+
ypos = ypos - panelheight;
795
+
entries(ii) = iUpdateEntry( entries(ii), true );
796
+
end
797
+
setappdata( fig, 'ProgressEntries', entries );
798
+
end % iRedraw
799
+
800
+
function cdata = iMakeCross( sz )
801
+
% Create a cross-shape icon of size sz*sz*3
802
+
803
+
cdata = diag(ones(sz,1),0) + diag(ones(sz-1,1),1) + diag(ones(sz-1,1),-1);
804
+
cdata = cdata + flip(cdata,2);
805
+
806
+
% Convert zeros to nans (transparent) and non-zeros to zero (black)
807
+
cdata(cdata == 0) = nan;
808
+
cdata(~isnan(cdata)) = 0;
809
+
810
+
% Convert to RGB
811
+
cdata = cat( 3, cdata, cdata, cdata );
812
+
end % iMakeCross
813
+
814
+
815
+
function iTimerFcn(fig)
816
+
% Timer callback for updating stuff every so often
817
+
entries = getappdata( fig, 'ProgressEntries' );
818
+
for ii=1:numel(entries)
819
+
if entries(ii).Busy
820
+
entries(ii) = iUpdateBusyEntry(entries(ii));
821
+
end
822
+
end
823
+
setappdata( fig, 'ProgressEntries', entries );
824
+
end % iTimerFcn
+21
include/representaAbsorcion.m
+21
include/representaAbsorcion.m
···
1
+
function representaAbsorcion(handles)
2
+
3
+
axes(handles.plotfig)
4
+
Absor = handles.SPLm{handles.Index,11};
5
+
Frecuencia = [125,250,500,1000,2000,4000,8000];
6
+
Fuente = string(handles.SPLm(handles.Index,2));
7
+
% Pasar de tercios de octava a octavas
8
+
Octavas = zeros(1,7);
9
+
a = 0;
10
+
for i=2:3:21
11
+
a = a+1;
12
+
Octavas(a) = (Absor(i-1)+Absor(i)+Absor(i+1))/3;
13
+
end
14
+
15
+
bar(Octavas*100)
16
+
set(gca,'XTickLabel', Frecuencia);
17
+
18
+
title('Coeficientes de absorci�n medio')
19
+
xlabel('Frecuencia (Hz)')
20
+
ylabel('\alpha (%)')
21
+
end
+21
include/representaDirectividad.m
+21
include/representaDirectividad.m
···
1
+
function representaDirectividad(handles)
2
+
3
+
axes(handles.plotfig)
4
+
Resp = handles.SPLm{handles.Index,14};
5
+
Frecuencia = [125,250,500,1000,2000,4000,8000];
6
+
Fuente = string(handles.SPLm(handles.Index,2));
7
+
% Obtener solo el valor de las octavas
8
+
Octavas = zeros(1,7);
9
+
a = 0;
10
+
for i=2:3:21
11
+
a = a+1;
12
+
Octavas(a) = Resp(i);
13
+
end
14
+
15
+
bar(Octavas)
16
+
set(gca,'XTickLabel', Frecuencia);
17
+
18
+
title(sprintf('Respuesta en frecuencia de la fuente %s',Fuente))
19
+
xlabel('Frecuencia (Hz)')
20
+
ylabel('DI (dB)')
21
+
end
+21
include/representaRT.m
+21
include/representaRT.m
···
1
+
function representaRT(handles)
2
+
3
+
axes(handles.plotfig)
4
+
RT = handles.SPLm{handles.Index,10};
5
+
Frecuencia = [125,250,500,1000,2000,4000,8000];
6
+
Fuente = string(handles.SPLm(handles.Index,2));
7
+
% Pasar de tercios de octava a octavas
8
+
Octavas = zeros(1,7);
9
+
a = 0;
10
+
for i=2:3:21
11
+
a = a+1;
12
+
Octavas(a) = (RT(i-1)+RT(i)+RT(i+1))/3;
13
+
end
14
+
15
+
bar(Octavas)
16
+
set(gca,'XTickLabel', Frecuencia);
17
+
18
+
title(sprintf('Tiempo de reverberaci�n con la fuente %s',Fuente))
19
+
xlabel('Frecuencia (Hz)')
20
+
ylabel('Tiempo de reverberaci�n (s)')
21
+
end
+21
include/representaRespuesta.m
+21
include/representaRespuesta.m
···
1
+
function representaRespuesta(handles)
2
+
3
+
axes(handles.plotfig)
4
+
Resp = handles.SPLm{handles.Index,13};
5
+
Frecuencia = [125,250,500,1000,2000,4000,8000];
6
+
Fuente = string(handles.SPLm(handles.Index,2));
7
+
% Pasar de tercios de octava a octavas
8
+
Octavas = zeros(1,7);
9
+
a = 0;
10
+
for i=2:3:21
11
+
a = a+1;
12
+
Octavas(a) = 10*log10(sum(10.^(Resp(i-1)/10)+10.^(Resp(i)/10)+10.^(Resp(i+1)/10)));
13
+
end
14
+
15
+
bar(Octavas)
16
+
set(gca,'XTickLabel', Frecuencia);
17
+
18
+
title(sprintf('Respuesta en frecuencia de la fuente %s',Fuente))
19
+
xlabel('Frecuencia (Hz)')
20
+
ylabel('Nivel SPL (dB)')
21
+
end
+331
include/representarCt.m
+331
include/representarCt.m
···
1
+
function representarCt(hObject,handles)
2
+
% Reprsentacion del parametro claridad
3
+
axes(handles.plotfig)
4
+
5
+
Rango = handles.Rango;
6
+
SPLmTot = handles.SPLm(:,19); % Matriz con los rayos sin procesar
7
+
SPLm = handles.SPLm(:,20); % Matriz con valores cada ms y por octavas
8
+
Distancia = cell2mat(handles.SPLm(:,6))'; % Vector de distancias fuente-receptor
9
+
Distplot = min(Distancia):0.01:max(Distancia); % Vector de distancias para las curvas
10
+
Fuente = string(handles.SPLm(:,2)');
11
+
Receptor = string(handles.SPLm(:,3)');
12
+
13
+
14
+
% Comprobar que bandas de frecuencias se han elegido
15
+
Bandas = [];
16
+
fsel = [];
17
+
if handles.frec125.Value==1
18
+
Bandas = [Bandas,1];
19
+
fsel = [fsel,125];
20
+
end
21
+
if handles.frec250.Value==1
22
+
Bandas = [Bandas,2];
23
+
fsel = [fsel,250];
24
+
end
25
+
if handles.frec500.Value==1
26
+
Bandas = [Bandas,3];
27
+
fsel = [fsel,500];
28
+
end
29
+
if handles.frec1000.Value==1
30
+
Bandas = [Bandas,4];
31
+
fsel = [fsel,1000];
32
+
end
33
+
if handles.frec2000.Value==1
34
+
Bandas = [Bandas,5];
35
+
fsel = [fsel,2000];
36
+
end
37
+
if handles.frec4000.Value==1
38
+
Bandas = [Bandas,6];
39
+
fsel = [fsel,4000];
40
+
end
41
+
if handles.frec8000.Value==1
42
+
Bandas = [Bandas,7];
43
+
fsel = [fsel,8000];
44
+
end
45
+
if isempty(Bandas); Bandas = 1:7; end
46
+
if isempty(fsel); fsel = 10^3*(2.^((-3:3))); end
47
+
48
+
%% Procesado de datos
49
+
% Separacion e integracion de los niveles antes y despues del rango
50
+
Dneto = zeros(1,numel(SPLm));
51
+
Eneto = zeros(1,numel(SPLm));
52
+
Lneto = zeros(1,numel(SPLm));
53
+
for i=1:numel(SPLm)
54
+
% Direct - Early - Late EXTRACTION
55
+
Dbruto = SPLm{i}(1,Bandas); % Niveles de 0 a 1 ms
56
+
Ebruto = SPLm{i}(2:Rango,Bandas); % Niveles de 1 ms a tiempo de integracion
57
+
Lbruto = SPLm{i}(Rango+1:end,Bandas); % Niveles de tiempo de integracion a infinito
58
+
59
+
% Ponderacion A
60
+
% f = fsel; % Frecuencias seleccionadas
61
+
% ponderacionA = 20*log10(...
62
+
% 12194^2*f.^4 ./...
63
+
% ((f.^2+20.6^2).*sqrt((f.^2+107.7^2).*(f.^2+737.9^2)).*(f.^2+12194^2)))+2;
64
+
% Dbruto = Dbruto + ponderacionA;
65
+
% Ebruto = Ebruto + ponderacionA;
66
+
% Lbruto = Lbruto + ponderacionA;
67
+
68
+
% Suma de niveles para los dos rangos
69
+
Dneto(i) = 10*log10(sum(10.^(Dbruto(:)/10)));
70
+
Eneto(i) = 10*log10(sum(10.^(Ebruto(:)/10)));
71
+
Lneto(i) = 10*log10(sum(10.^(Lbruto(:)/10)));
72
+
end
73
+
74
+
% Campo util es la suma del campo directo y el early
75
+
Uneto = 10*log10(10.^(Dneto/10)+10.^(Eneto/10));
76
+
77
+
%% Representacion grafica del campo Util y el Perjudicial
78
+
79
+
%% Puntos y curva de 0ms a Valor
80
+
% Ordenar valores de distancia de menor a mayor manteniendo su valor de SPl
81
+
% asociado
82
+
Dist_U = sortrows([Distancia',Uneto']);
83
+
% Si tiene distancias duplicadas/iguales se promedia el nivel de las dos
84
+
% distancias iguales y elimina el duplicado
85
+
[C,~,idx] = unique(Dist_U(:,1),'stable');
86
+
val = accumarray(idx,Dist_U(:,2),[],@(x) 10*log10(mean(10.^(x/10))));
87
+
Dist_U = [C,val];
88
+
% A�adir por extrapolacion valor en la posicion de 0 metros hasta la posicion del
89
+
% receptor mas cercano
90
+
%V0_0toVal=interp1(Dist_U(:,1),Dist_U(:,2),0,'linear','extrap');
91
+
%Dist_U = [[0;V0_0toVal]';Dist_U];
92
+
% Separar los valores para graficar y convertirlo en vector fila
93
+
Dist = Dist_U(:,1);
94
+
Urec = Dist_U(:,2);
95
+
% Crear curva potencial. Obtiene coeficientes (Ucoef) y datos de ajuste
96
+
% (Ures)
97
+
[Ucoef,Ures,Uoutput]=fit(Dist,Urec,'power1');
98
+
Uplot=Ucoef.a*Distplot.^Ucoef.b;
99
+
100
+
101
+
102
+
%% Puntos y curva de valor a infinito ms
103
+
% Ordenar valores de distancia de menor a mayor manteniendo su valor de SPl
104
+
% asociado
105
+
Dist_P = sortrows([Distancia',Lneto']);
106
+
% Si tiene distancias duplicadas/iguales se promedia el nivel de las dos
107
+
% distancias iguales y elimina el duplicado
108
+
[C,~,idx] = unique(Dist_P(:,1),'stable');
109
+
val = accumarray(idx,Dist_P(:,2),[],@(x) 10*log10(mean(10.^(x/10))));
110
+
Dist_P = [C,val];
111
+
% A�adir por extrapolacion valor en la posicion de 0 metros hasta la posicion del
112
+
% receptor mas cercano
113
+
%V0_0toVal=interp1(Dist_P(:,1),Dist_P(:,2),0,'linear','extrap');
114
+
%Dist_P = [[0;V0_0toVal]';Dist_P];
115
+
% Separar los valores para graficar y convertirlo en vector fila
116
+
Dist = Dist_P(:,1);
117
+
Prec = Dist_P(:,2);
118
+
% Obtener coeificentes curva polinomica de grado 1
119
+
[Pcoef,Pres,Poutput] = fit(Dist,Prec,'poly1');
120
+
% Crea la curva
121
+
Pplot = Pcoef.p1*Distplot+Pcoef.p2;
122
+
123
+
124
+
% Claridad EASE
125
+
% Muestra los puntos
126
+
plot(Dist,Urec-Prec,'^b','MarkerFaceColor','b'),hold on
127
+
% Muestra la curva
128
+
Curv(1)=plot(Distplot,Uplot-Pplot,'b','LineWidth',1.2);
129
+
130
+
% A�ade las rejillas de las graficas
131
+
grid on
132
+
grid minor
133
+
134
+
% Leyenda
135
+
leyenda{1} = sprintf('C_{%s} - EASE',num2str(Rango));
136
+
137
+
%% Teoria corregida
138
+
% Si se ha marcado el calculo de la teoria revisada corregida se calcula y
139
+
% se a�ade a la grafica
140
+
if get(handles.teoriacorregida,'value')==1
141
+
%%% PARAMETROS
142
+
% Nivel de presion total de la fuente a 1 metro
143
+
SPL = 10*log10(sum(10.^(handles.SPLm{handles.Index,13}/10)));
144
+
% Volumen (m3)
145
+
V = handles.SPLm{handles.Index,15};
146
+
% Superficie (m2)
147
+
S = handles.SPLm{handles.Index,16};
148
+
% Mean Free Path
149
+
mfp = 4*V/S;
150
+
% Absorcion media
151
+
alpha = mean(handles.SPLm{handles.Index,11});
152
+
% Absorci�n equivalente de Eyring
153
+
A = -S*log(1-alpha);
154
+
% Temperatura de la sala
155
+
temp = handles.SPLm{handles.Index,17};
156
+
% Velocidad del sonido
157
+
c = 331.4 + 0.6*temp;
158
+
% Densidad del aire. Presion Atm / (constante gas*Temp
159
+
% absoluta)
160
+
rho = 101325 / (287.058*(273.15+temp));
161
+
% Impedancia acustica del aire
162
+
Z = rho*c;
163
+
% Presion de referencia
164
+
pref = 2*10^(-5);
165
+
% Factor de conversion
166
+
Fac = Z/(pref^2);
167
+
% Directividad de la fuente
168
+
Q = str2double(get(handles.Qdir,'string'));
169
+
% Nivel de potencia de la fuente (SPL a 1 metro). Se supone Q=1
170
+
SWL = SPL - 10*log10(1/(4*pi*1));
171
+
% Potencia de la fuente en vatios
172
+
W = 10.^(SWL/10)*10^-12;
173
+
% Tiempo de reverberacion de Eyring
174
+
T = 6*log(10)/c * mfp * 1/(-log(1-alpha));
175
+
% Rango de integracion temporal (s)
176
+
t = handles.Rango * 10^(-3);
177
+
178
+
%% Obtencion de la curvas directo y early
179
+
%% Curva del campo directo
180
+
% Ordenar valores de distancia de menor a mayor manteniendo su valor de SPl
181
+
% asociado
182
+
Dist_D = sortrows([Distancia',Dneto']);
183
+
% Si tiene distancias duplicadas/iguales se promedia el nivel de las dos
184
+
% distancias iguales y elimina el duplicado
185
+
[C,~,idx] = unique(Dist_D(:,1),'stable');
186
+
val = accumarray(idx,Dist_D(:,2),[],@(x) 10*log10(mean(10.^(x/10))));
187
+
Dist_D = [C,val];
188
+
% A�adir por extrapolacion valor en la posicion de 0 metros hasta la posicion del
189
+
% receptor mas cercano
190
+
%V0_0toVal=interp1(Dist_U(:,1),Dist_U(:,2),0,'linear','extrap');
191
+
%Dist_U = [[0;V0_0toVal]';Dist_U];
192
+
% Separar los valores para graficar y convertirlo en vector fila
193
+
Dist = Dist_D(:,1);
194
+
Drec = Dist_D(:,2);
195
+
% Crear curva potencial. Obtiene coeficientes (Ucoef) y datos de ajuste
196
+
% (Ures)
197
+
[Dcoef,Dres]=fit(Dist,Drec,'power1');
198
+
Dplot=Dcoef.a*Distplot.^Dcoef.b;
199
+
200
+
%% Curva del campo early
201
+
% Ordenar valores de distancia de menor a mayor manteniendo su valor de SPl
202
+
% asociado
203
+
Dist_E = sortrows([Distancia',Eneto']);
204
+
% Si tiene distancias duplicadas/iguales se promedia el nivel de las dos
205
+
% distancias iguales y elimina el duplicado
206
+
[C,~,idx] = unique(Dist_E(:,1),'stable');
207
+
val = accumarray(idx,Dist_E(:,2),[],@(x) 10*log10(mean(10.^(x/10))));
208
+
Dist_E = [C,val];
209
+
% A�adir por extrapolacion valor en la posicion de 0 metros hasta la posicion del
210
+
% receptor mas cercano
211
+
%V0_0toVal=interp1(Dist_U(:,1),Dist_U(:,2),0,'linear','extrap');
212
+
%Dist_U = [[0;V0_0toVal]';Dist_U];
213
+
% Separar los valores para graficar y convertirlo en vector fila
214
+
Dist = Dist_E(:,1);
215
+
Erec = Dist_E(:,2);
216
+
% Crear curva potencial. Obtiene coeficientes (Ucoef) y datos de ajuste
217
+
% (Ures)
218
+
[Ecoef,Eres]=fit(Dist,Erec,'power1');
219
+
Eplot=Ecoef.a*Distplot.^Ecoef.b;
220
+
221
+
% Para no mostrar errores cuando se repite el calculo de coeficientes
222
+
warning ('off','all');
223
+
224
+
%% Calculo de coeficiente para el campo directo
225
+
for ii = 1:str2double(get(handles.numIntentos,'string')) % 100 intentos de obtener los coeficientes
226
+
try % Intenta el calculo con puntos de inicio random
227
+
opD = fitoptions('Method','NonlinearLeastSquares',...
228
+
'Robust','LAR',...
229
+
'Lower',[0,0],... % Minimo valor de los coeficientes
230
+
'Upper',[20,20],... % Maximo valor de los coeficientes
231
+
'MaxFunEvals',3000,...
232
+
'MaxIter',3000,...
233
+
'TolFun',10^-3);
234
+
fD = fittype('Fac * W*Q./(4*pi*dist.^2)*Cd',...
235
+
'dependent',{'y'},...
236
+
'independent',{'dist'},...
237
+
'problem', {'Q','W','Fac'},...
238
+
'coefficients',{'Cd'},...
239
+
'options',opD);
240
+
[fitiD,gofD,outputD] = fit(Distplot',10.^(Dplot'/10),fD,'problem',{Q,W,Fac});
241
+
break;
242
+
catch % Si no se ha podido obtener unos coeficientes ejecuta lo siguiente
243
+
if ii==str2double(get(handles.numIntentos,'string'))
244
+
hold off
245
+
ed = errordlg('No se ha podido obtener coeficientes (Directo), intentalo de nuevo','Error');
246
+
set(ed, 'WindowStyle', 'modal');
247
+
uiwait(ed);
248
+
return;
249
+
end
250
+
end
251
+
end
252
+
253
+
254
+
%% Calculo de coeficientes epsilon y Cl para campo late
255
+
for ii = 1:str2double(get(handles.numIntentos,'string')) % 100 intentos de obtener los coeficientes
256
+
try % Intenta el calculo con puntos de inicio random
257
+
fL = fittype('Fac * (4*W)/(S*(-log(1-alpha))) * exp(-(13.82*(dist/c+t)*el/T)) * Cl',...
258
+
'dependent',{'y'},...
259
+
'independent',{'dist'},...
260
+
'problem', {'W','S','alpha','t','T','c','Fac'},...
261
+
'coefficients',{'el','Cl'});
262
+
[fitiL,gofL,outputL] = fit(Distplot',10.^(Pplot'/10),fL,'problem',{W,S,alpha,t,T,c,Fac});
263
+
break;
264
+
catch % Si no se ha podido obtener unos coeficientes ejecuta lo siguiente
265
+
if ii==str2double(get(handles.numIntentos,'string'))
266
+
hold off
267
+
ed = errordlg('No se ha podido obtener coeficientes (Late), intentalo de nuevo','Error');
268
+
set(ed, 'WindowStyle', 'modal');
269
+
uiwait(ed);
270
+
return;
271
+
end
272
+
end
273
+
end
274
+
275
+
%% Calculo de coeficientes epsilon y Ce para campo early
276
+
for ii = 1:str2double(get(handles.numIntentos,'string')) % 100 intentos de obtener los coeficientes
277
+
try % Intenta el calculo con puntos de inicio random
278
+
opE = fitoptions('Method','NonlinearLeastSquares',...
279
+
'Robust','LAR',...
280
+
'MaxFunEvals',3000,...
281
+
'MaxIter',3000,...
282
+
'TolFun',10^-3);
283
+
fE = fittype('Fac * ((4*W)./(S*(-log(1-alpha))*dist) .* (exp(-(13.82*(dist/c)*ee/T))*Ce - exp(-(13.82*(dist/c+t)*el/T)) * Cl))',...
284
+
'dependent',{'y'},...
285
+
'independent',{'dist'},...
286
+
'problem', {'W','S','alpha','T','t','c','Fac','el','Cl'},...
287
+
'coefficients',{'ee','Ce'},...
288
+
'options',opE);
289
+
[fitiE,gofE,outputE] = fit(Distplot',10.^(Eplot'/10),fE,'problem',{W,S,alpha,T,t,c,Fac,fitiL.el,fitiL.Cl});
290
+
break;
291
+
catch % Si no se ha podido obtener unos coeficientes ejecuta lo siguiente
292
+
if ii==str2double(get(handles.numIntentos,'string'))
293
+
hold off
294
+
ed = errordlg('No se ha podido obtener coeficientes (Early), intentalo de nuevo','Error');
295
+
set(ed, 'WindowStyle', 'modal');
296
+
uiwait(ed);
297
+
return;
298
+
end
299
+
end
300
+
end
301
+
warning ('on','all');
302
+
303
+
PDplot = 10*log10(Fac * ( W*Q./(4*pi*Distplot.^2)*fitiD.Cd));
304
+
Peplot = 10*log10(Fac *(...
305
+
(4*W)./(S*(-log(1-alpha))*Distplot) .* ...
306
+
(exp(-(13.82*(Distplot/c)*fitiE.ee/T))*fitiE.Ce - exp(-(13.82*(Distplot/c+t)*fitiL.el/T))*fitiL.Cl)));
307
+
Pdeplot = 10*log10(10.^(PDplot/10)+10.^(Peplot/10));
308
+
Plateplot = 10*log10(Fac .* (4*W)./(S*(-log(1-alpha))) .* exp(-(13.82.*(Distplot./c+t).*fitiL.el/T)) .* fitiL.Cl);
309
+
310
+
Curv(2) = plot(Distplot,Pdeplot - Plateplot,'-.r','LineWidth',1);
311
+
312
+
% Leyendas
313
+
leyenda{2} = sprintf('C_{%s} - Te�rico',num2str(Rango));
314
+
315
+
end
316
+
317
+
lgdw=legend(Curv,leyenda);
318
+
lgdw.FontSize = 11;
319
+
%lgdw.Location = 'northeastoutside';
320
+
% Si se utiliza Matlab 2014b o anterior no a�ade el titulo a la leyenda ya
321
+
% que no es compatible
322
+
if ~verLessThan('matlab','8.7')
323
+
title(lgdw,'Curvas')
324
+
end
325
+
ylabel(sprintf('C_{%s} (dB)',num2str(Rango)))
326
+
xlabel('Distancia (m)')
327
+
title({sprintf('C_{%s}',num2str(Rango));...
328
+
['\color{blue} ','Fuente: ',sprintf(': %s',Fuente(1))]});
329
+
330
+
hold off
331
+
+330
include/representarDt.m
+330
include/representarDt.m
···
1
+
function representarDt(hObject,handles)
2
+
% Reprsentacion del parametro definicion
3
+
axes(handles.plotfig)
4
+
Rango = handles.Rango;
5
+
SPLmTot = handles.SPLm(:,19); % Matriz con los rayos sin procesar
6
+
SPLm = handles.SPLm(:,20); % Matriz con valores cada ms y por octavas
7
+
Distancia = cell2mat(handles.SPLm(:,6))'; % Vector de distancias fuente-receptor
8
+
Distplot = min(Distancia):0.01:max(Distancia); % Vector de distancias para las curvas
9
+
Fuente = string(handles.SPLm(:,2)');
10
+
Receptor = string(handles.SPLm(:,3)');
11
+
12
+
13
+
% Comprobar que bandas de frecuencias se han elegido
14
+
Bandas = [];
15
+
fsel = [];
16
+
if handles.frec125.Value==1
17
+
Bandas = [Bandas,1];
18
+
fsel = [fsel,125];
19
+
end
20
+
if handles.frec250.Value==1
21
+
Bandas = [Bandas,2];
22
+
fsel = [fsel,250];
23
+
end
24
+
if handles.frec500.Value==1
25
+
Bandas = [Bandas,3];
26
+
fsel = [fsel,500];
27
+
end
28
+
if handles.frec1000.Value==1
29
+
Bandas = [Bandas,4];
30
+
fsel = [fsel,1000];
31
+
end
32
+
if handles.frec2000.Value==1
33
+
Bandas = [Bandas,5];
34
+
fsel = [fsel,2000];
35
+
end
36
+
if handles.frec4000.Value==1
37
+
Bandas = [Bandas,6];
38
+
fsel = [fsel,4000];
39
+
end
40
+
if handles.frec8000.Value==1
41
+
Bandas = [Bandas,7];
42
+
fsel = [fsel,8000];
43
+
end
44
+
if isempty(Bandas); Bandas = 1:7; end
45
+
if isempty(fsel); fsel = 10^3*(2.^((-3:3))); end
46
+
47
+
%% Procesado de datos
48
+
% Separacion e integracion de los niveles antes y despues del rango
49
+
Dneto = zeros(1,numel(SPLm));
50
+
Eneto = zeros(1,numel(SPLm));
51
+
Lneto = zeros(1,numel(SPLm));
52
+
for i=1:numel(SPLm)
53
+
% Direct - Early - Late EXTRACTION
54
+
Dbruto = SPLm{i}(1,Bandas); % Niveles de 0 a 1 ms
55
+
Ebruto = SPLm{i}(2:Rango,Bandas); % Niveles de 1 ms a tiempo de integracion
56
+
Lbruto = SPLm{i}(Rango+1:end,Bandas); % Niveles de tiempo de integracion a infinito
57
+
58
+
% Ponderacion A
59
+
% f = fsel; % Frecuencias seleccionadas
60
+
% ponderacionA = 20*log10(...
61
+
% 12194^2*f.^4 ./...
62
+
% ((f.^2+20.6^2).*sqrt((f.^2+107.7^2).*(f.^2+737.9^2)).*(f.^2+12194^2)))+2;
63
+
% Dbruto = Dbruto + ponderacionA;
64
+
% Ebruto = Ebruto + ponderacionA;
65
+
% Lbruto = Lbruto + ponderacionA;
66
+
67
+
% Suma de niveles para los dos rangos
68
+
Dneto(i) = 10*log10(sum(10.^(Dbruto(:)/10)));
69
+
Eneto(i) = 10*log10(sum(10.^(Ebruto(:)/10)));
70
+
Lneto(i) = 10*log10(sum(10.^(Lbruto(:)/10)));
71
+
end
72
+
73
+
% Campo util es la suma del campo directo y el early
74
+
Uneto = 10*log10(10.^(Dneto/10)+10.^(Eneto/10));
75
+
76
+
%% Representacion grafica del campo Util y el Perjudicial
77
+
78
+
%% Puntos y curva de 0ms a Valor
79
+
% Ordenar valores de distancia de menor a mayor manteniendo su valor de SPl
80
+
% asociado
81
+
Dist_U = sortrows([Distancia',Uneto']);
82
+
% Si tiene distancias duplicadas/iguales se promedia el nivel de las dos
83
+
% distancias iguales y elimina el duplicado
84
+
[C,~,idx] = unique(Dist_U(:,1),'stable');
85
+
val = accumarray(idx,Dist_U(:,2),[],@(x) 10*log10(mean(10.^(x/10))));
86
+
Dist_U = [C,val];
87
+
% A�adir por extrapolacion valor en la posicion de 0 metros hasta la posicion del
88
+
% receptor mas cercano
89
+
%V0_0toVal=interp1(Dist_U(:,1),Dist_U(:,2),0,'linear','extrap');
90
+
%Dist_U = [[0;V0_0toVal]';Dist_U];
91
+
% Separar los valores para graficar y convertirlo en vector fila
92
+
Dist = Dist_U(:,1);
93
+
Urec = Dist_U(:,2);
94
+
% Crear curva potencial. Obtiene coeficientes (Ucoef) y datos de ajuste
95
+
% (Ures)
96
+
[Ucoef,Ures,Uoutput]=fit(Dist,Urec,'power1');
97
+
Uplot=Ucoef.a*Distplot.^Ucoef.b;
98
+
99
+
100
+
101
+
%% Puntos y curva de valor a infinito ms
102
+
% Ordenar valores de distancia de menor a mayor manteniendo su valor de SPl
103
+
% asociado
104
+
Dist_P = sortrows([Distancia',Lneto']);
105
+
% Si tiene distancias duplicadas/iguales se promedia el nivel de las dos
106
+
% distancias iguales y elimina el duplicado
107
+
[C,~,idx] = unique(Dist_P(:,1),'stable');
108
+
val = accumarray(idx,Dist_P(:,2),[],@(x) 10*log10(mean(10.^(x/10))));
109
+
Dist_P = [C,val];
110
+
% A�adir por extrapolacion valor en la posicion de 0 metros hasta la posicion del
111
+
% receptor mas cercano
112
+
%V0_0toVal=interp1(Dist_P(:,1),Dist_P(:,2),0,'linear','extrap');
113
+
%Dist_P = [[0;V0_0toVal]';Dist_P];
114
+
% Separar los valores para graficar y convertirlo en vector fila
115
+
Dist = Dist_P(:,1);
116
+
Prec = Dist_P(:,2);
117
+
% Obtener coeificentes curva polinomica de grado 1
118
+
[Pcoef,Pres,Poutput] = fit(Dist,Prec,'poly1');
119
+
% Crea la curva
120
+
Pplot = Pcoef.p1*Distplot+Pcoef.p2;
121
+
122
+
123
+
% Claridad EASE
124
+
% Muestra los puntos
125
+
plot(Dist,10.^(Urec/10)./( 10.^(Urec/10) + 10.^(Prec/10) ),'^b','MarkerFaceColor','b'),hold on
126
+
% Muestra la curva
127
+
Curv(1)=plot(Distplot,10.^(Uplot/10)./( 10.^(Uplot/10) + 10.^(Pplot/10) ),'b','LineWidth',1.2);
128
+
129
+
% A�ade las rejillas de las graficas
130
+
grid on
131
+
grid minor
132
+
133
+
% Leyenda
134
+
leyenda{1} = sprintf('D_{%s} - EASE',num2str(Rango));
135
+
136
+
%% Teoria corregida
137
+
% Si se ha marcado el calculo de la teoria revisada corregida se calcula y
138
+
% se a�ade a la grafica
139
+
if get(handles.teoriacorregida,'value')==1
140
+
%%% PARAMETROS
141
+
% Nivel de presion total de la fuente a 1 metro
142
+
SPL = 10*log10(sum(10.^(handles.SPLm{handles.Index,13}/10)));
143
+
% Volumen (m3)
144
+
V = handles.SPLm{handles.Index,15};
145
+
% Superficie (m2)
146
+
S = handles.SPLm{handles.Index,16};
147
+
% Mean Free Path
148
+
mfp = 4*V/S;
149
+
% Absorcion media
150
+
alpha = mean(handles.SPLm{handles.Index,11});
151
+
% Absorci�n equivalente de Eyring
152
+
A = -S*log(1-alpha);
153
+
% Temperatura de la sala
154
+
temp = handles.SPLm{handles.Index,17};
155
+
% Velocidad del sonido
156
+
c = 331.4 + 0.6*temp;
157
+
% Densidad del aire. Presion Atm / (constante gas*Temp
158
+
% absoluta)
159
+
rho = 101325 / (287.058*(273.15+temp));
160
+
% Impedancia acustica del aire
161
+
Z = rho*c;
162
+
% Presion de referencia
163
+
pref = 2*10^(-5);
164
+
% Factor de conversion
165
+
Fac = Z/(pref^2);
166
+
% Directividad de la fuente
167
+
Q = str2double(get(handles.Qdir,'string'));
168
+
% Nivel de potencia de la fuente (SPL a 1 metro). Se supone Q=1
169
+
SWL = SPL - 10*log10(1/(4*pi*1));
170
+
% Potencia de la fuente en vatios
171
+
W = 10.^(SWL/10)*10^-12;
172
+
% Tiempo de reverberacion de Eyring
173
+
T = 6*log(10)/c * mfp * 1/(-log(1-alpha));
174
+
% Rango de integracion temporal (s)
175
+
t = handles.Rango * 10^(-3);
176
+
177
+
%% Obtencion de la curvas directo y early
178
+
%% Curva del campo directo
179
+
% Ordenar valores de distancia de menor a mayor manteniendo su valor de SPl
180
+
% asociado
181
+
Dist_D = sortrows([Distancia',Dneto']);
182
+
% Si tiene distancias duplicadas/iguales se promedia el nivel de las dos
183
+
% distancias iguales y elimina el duplicado
184
+
[C,~,idx] = unique(Dist_D(:,1),'stable');
185
+
val = accumarray(idx,Dist_D(:,2),[],@(x) 10*log10(mean(10.^(x/10))));
186
+
Dist_D = [C,val];
187
+
% A�adir por extrapolacion valor en la posicion de 0 metros hasta la posicion del
188
+
% receptor mas cercano
189
+
%V0_0toVal=interp1(Dist_U(:,1),Dist_U(:,2),0,'linear','extrap');
190
+
%Dist_U = [[0;V0_0toVal]';Dist_U];
191
+
% Separar los valores para graficar y convertirlo en vector fila
192
+
Dist = Dist_D(:,1);
193
+
Drec = Dist_D(:,2);
194
+
% Crear curva potencial. Obtiene coeficientes (Ucoef) y datos de ajuste
195
+
% (Ures)
196
+
[Dcoef,Dres]=fit(Dist,Drec,'power1');
197
+
Dplot=Dcoef.a*Distplot.^Dcoef.b;
198
+
199
+
%% Curva del campo early
200
+
% Ordenar valores de distancia de menor a mayor manteniendo su valor de SPl
201
+
% asociado
202
+
Dist_E = sortrows([Distancia',Eneto']);
203
+
% Si tiene distancias duplicadas/iguales se promedia el nivel de las dos
204
+
% distancias iguales y elimina el duplicado
205
+
[C,~,idx] = unique(Dist_E(:,1),'stable');
206
+
val = accumarray(idx,Dist_E(:,2),[],@(x) 10*log10(mean(10.^(x/10))));
207
+
Dist_E = [C,val];
208
+
% A�adir por extrapolacion valor en la posicion de 0 metros hasta la posicion del
209
+
% receptor mas cercano
210
+
%V0_0toVal=interp1(Dist_U(:,1),Dist_U(:,2),0,'linear','extrap');
211
+
%Dist_U = [[0;V0_0toVal]';Dist_U];
212
+
% Separar los valores para graficar y convertirlo en vector fila
213
+
Dist = Dist_E(:,1);
214
+
Erec = Dist_E(:,2);
215
+
% Crear curva potencial. Obtiene coeficientes (Ucoef) y datos de ajuste
216
+
% (Ures)
217
+
[Ecoef,Eres]=fit(Dist,Erec,'power1');
218
+
Eplot=Ecoef.a*Distplot.^Ecoef.b;
219
+
220
+
% Para no mostrar errores cuando se repite el calculo de coeficientes
221
+
warning ('off','all');
222
+
223
+
%% Calculo de coeficiente para el campo directo
224
+
for ii = 1:str2double(get(handles.numIntentos,'string')) % 100 intentos de obtener los coeficientes
225
+
try % Intenta el calculo con puntos de inicio random
226
+
opD = fitoptions('Method','NonlinearLeastSquares',...
227
+
'Robust','LAR',...
228
+
'Lower',[0,0],... % Minimo valor de los coeficientes
229
+
'Upper',[20,20],... % Maximo valor de los coeficientes
230
+
'MaxFunEvals',3000,...
231
+
'MaxIter',3000,...
232
+
'TolFun',10^-3);
233
+
fD = fittype('Fac * W*Q./(4*pi*dist.^2)*Cd',...
234
+
'dependent',{'y'},...
235
+
'independent',{'dist'},...
236
+
'problem', {'Q','W','Fac'},...
237
+
'coefficients',{'Cd'},...
238
+
'options',opD);
239
+
[fitiD,gofD,outputD] = fit(Distplot',10.^(Dplot'/10),fD,'problem',{Q,W,Fac});
240
+
break;
241
+
catch % Si no se ha podido obtener unos coeficientes ejecuta lo siguiente
242
+
if ii==str2double(get(handles.numIntentos,'string'))
243
+
hold off
244
+
ed = errordlg('No se ha podido obtener coeficientes (Directo), intentalo de nuevo','Error');
245
+
set(ed, 'WindowStyle', 'modal');
246
+
uiwait(ed);
247
+
return;
248
+
end
249
+
end
250
+
end
251
+
252
+
253
+
%% Calculo de coeficientes epsilon y Cl para campo late
254
+
for ii = 1:str2double(get(handles.numIntentos,'string')) % 100 intentos de obtener los coeficientes
255
+
try % Intenta el calculo con puntos de inicio random
256
+
fL = fittype('Fac * (4*W)/(S*(-log(1-alpha))) * exp(-(13.82*(dist/c+t)*el/T)) * Cl',...
257
+
'dependent',{'y'},...
258
+
'independent',{'dist'},...
259
+
'problem', {'W','S','alpha','t','T','c','Fac'},...
260
+
'coefficients',{'el','Cl'});
261
+
[fitiL,gofL,outputL] = fit(Distplot',10.^(Pplot'/10),fL,'problem',{W,S,alpha,t,T,c,Fac});
262
+
break;
263
+
catch % Si no se ha podido obtener unos coeficientes ejecuta lo siguiente
264
+
if ii==str2double(get(handles.numIntentos,'string'))
265
+
hold off
266
+
ed = errordlg('No se ha podido obtener coeficientes (Late), intentalo de nuevo','Error');
267
+
set(ed, 'WindowStyle', 'modal');
268
+
uiwait(ed);
269
+
return;
270
+
end
271
+
end
272
+
end
273
+
274
+
%% Calculo de coeficientes epsilon y Ce para campo early
275
+
for ii = 1:str2double(get(handles.numIntentos,'string')) % 100 intentos de obtener los coeficientes
276
+
try % Intenta el calculo con puntos de inicio random
277
+
opE = fitoptions('Method','NonlinearLeastSquares',...
278
+
'Robust','LAR',...
279
+
'MaxFunEvals',3000,...
280
+
'MaxIter',3000,...
281
+
'TolFun',10^-3);
282
+
fE = fittype('Fac * ((4*W)./(S*(-log(1-alpha))*dist) .* (exp(-(13.82*(dist/c)*ee/T))*Ce - exp(-(13.82*(dist/c+t)*el/T)) * Cl))',...
283
+
'dependent',{'y'},...
284
+
'independent',{'dist'},...
285
+
'problem', {'W','S','alpha','T','t','c','Fac','el','Cl'},...
286
+
'coefficients',{'ee','Ce'},...
287
+
'options',opE);
288
+
[fitiE,gofE,outputE] = fit(Distplot',10.^(Eplot'/10),fE,'problem',{W,S,alpha,T,t,c,Fac,fitiL.el,fitiL.Cl});
289
+
break;
290
+
catch % Si no se ha podido obtener unos coeficientes ejecuta lo siguiente
291
+
if ii==str2double(get(handles.numIntentos,'string'))
292
+
hold off
293
+
ed = errordlg('No se ha podido obtener coeficientes (Early), intentalo de nuevo','Error');
294
+
set(ed, 'WindowStyle', 'modal');
295
+
uiwait(ed);
296
+
return;
297
+
end
298
+
end
299
+
end
300
+
warning ('on','all');
301
+
302
+
PDplot = 10*log10(Fac * ( W*Q./(4*pi*Distplot.^2)*fitiD.Cd));
303
+
Peplot = 10*log10(Fac *(...
304
+
(4*W)./(S*(-log(1-alpha))*Distplot) .* ...
305
+
(exp(-(13.82*(Distplot/c)*fitiE.ee/T))*fitiE.Ce - exp(-(13.82*(Distplot/c+t)*fitiL.el/T))*fitiL.Cl)));
306
+
Pdeplot = 10*log10(10.^(PDplot/10)+10.^(Peplot/10));
307
+
Plateplot = 10*log10(Fac .* (4*W)./(S*(-log(1-alpha))) .* exp(-(13.82.*(Distplot./c+t).*fitiL.el/T)) .* fitiL.Cl);
308
+
309
+
Curv(2) = plot(Distplot,10.^(Pdeplot/10)./( 10.^(Pdeplot/10) + 10.^(Plateplot/10) ),'-.r','LineWidth',1);
310
+
311
+
% Leyendas
312
+
leyenda{2} = sprintf('D_{%s} - Te�rico',num2str(Rango));
313
+
314
+
end
315
+
316
+
lgdw=legend(Curv,leyenda);
317
+
lgdw.FontSize = 11;
318
+
%lgdw.Location = 'northeastoutside';
319
+
% Si se utiliza Matlab 2014b o anterior no a�ade el titulo a la leyenda ya
320
+
% que no es compatible
321
+
if ~verLessThan('matlab','8.7')
322
+
title(lgdw,'Curvas')
323
+
end
324
+
ylabel(sprintf('D_{%s} (dB)',num2str(Rango)))
325
+
xlabel('Distancia (m)')
326
+
title({sprintf('D_{%s}',num2str(Rango));...
327
+
['\color{blue} ','Fuente: ',sprintf(': %s',Fuente(1))]});
328
+
329
+
hold off
330
+
+30
include/representarEcograma.m
+30
include/representarEcograma.m
···
1
+
function representarEcograma(handles)
2
+
3
+
axes(handles.plotfig)
4
+
Index = handles.Index;
5
+
SPLm = handles.SPLm(:,1);
6
+
SPLmTot = handles.SPLm(:,19);
7
+
Distancia = cell2mat(handles.SPLm(:,6))';
8
+
Fuente = string(handles.SPLm(:,2)');
9
+
Receptor = string(handles.SPLm(:,3)');
10
+
11
+
12
+
% Obtener ecograma del receptor Index, sumando las octavas
13
+
Ecogra = 10*log10(sum(10.^((SPLmTot{Index}(:,2:end)')/10)));
14
+
tiempo = SPLmTot{Index}(:,1);
15
+
16
+
17
+
18
+
%% Representacion grafica
19
+
20
+
21
+
% Ecograma
22
+
stem(tiempo,Ecogra,'Marker','none')
23
+
24
+
title({sprintf('Ecograma del receptor: %s',Receptor(Index));...
25
+
['\color{blue} ','Fuente: ',sprintf(': %s',Fuente(1))]})
26
+
xlabel('Tiempo (ms)')
27
+
ylabel('Nivel SPL (dB)')
28
+
% Limite del eje x con margen a ambos lados
29
+
xlim([-max(tiempo)/10,max(tiempo)+max(tiempo)/10])
30
+
end
+367
include/representarG.m
+367
include/representarG.m
···
1
+
function representarG(hObject,handles)
2
+
% Reprsentacion del parametro sonoridad
3
+
4
+
axes(handles.plotfig)
5
+
Rango = handles.Rango;
6
+
SPLmTot = handles.SPLm(:,19); % Matriz con los rayos sin procesar
7
+
SPLm = handles.SPLm(:,20); % Matriz con valores cada ms y por octavas
8
+
Distancia = cell2mat(handles.SPLm(:,6))'; % Vector de distancias fuente-receptor
9
+
Distplot = min(Distancia):0.01:max(Distancia); % Vector de distancias para las curvas
10
+
Fuente = string(handles.SPLm(:,2)');
11
+
Receptor = string(handles.SPLm(:,3)');
12
+
13
+
14
+
% Comprobar que bandas de frecuencias se han elegido
15
+
Bandas = [];
16
+
fsel = [];
17
+
if handles.frec125.Value==1
18
+
Bandas = [Bandas,1];
19
+
fsel = [fsel,125];
20
+
end
21
+
if handles.frec250.Value==1
22
+
Bandas = [Bandas,2];
23
+
fsel = [fsel,250];
24
+
end
25
+
if handles.frec500.Value==1
26
+
Bandas = [Bandas,3];
27
+
fsel = [fsel,500];
28
+
end
29
+
if handles.frec1000.Value==1
30
+
Bandas = [Bandas,4];
31
+
fsel = [fsel,1000];
32
+
end
33
+
if handles.frec2000.Value==1
34
+
Bandas = [Bandas,5];
35
+
fsel = [fsel,2000];
36
+
end
37
+
if handles.frec4000.Value==1
38
+
Bandas = [Bandas,6];
39
+
fsel = [fsel,4000];
40
+
end
41
+
if handles.frec8000.Value==1
42
+
Bandas = [Bandas,7];
43
+
fsel = [fsel,8000];
44
+
end
45
+
if isempty(Bandas); Bandas = 1:7; end
46
+
if isempty(fsel); fsel = 10^3*(2.^((-3:3))); end
47
+
48
+
%% Procesado de datos
49
+
% Separacion e integracion de los niveles antes y despues del rango
50
+
Dneto = zeros(1,numel(SPLm));
51
+
Eneto = zeros(1,numel(SPLm));
52
+
Lneto = zeros(1,numel(SPLm));
53
+
for i=1:numel(SPLm)
54
+
% Direct - Early - Late EXTRACTION
55
+
Dbruto = SPLm{i}(1,Bandas); % Niveles de 0 a 1 ms
56
+
Ebruto = SPLm{i}(2:Rango,Bandas); % Niveles de 1 ms a tiempo de integracion
57
+
Lbruto = SPLm{i}(Rango+1:end,Bandas); % Niveles de tiempo de integracion a infinito
58
+
59
+
% Ponderacion A
60
+
% f = fsel; % Frecuencias seleccionadas
61
+
% ponderacionA = 20*log10(...
62
+
% 12194^2*f.^4 ./...
63
+
% ((f.^2+20.6^2).*sqrt((f.^2+107.7^2).*(f.^2+737.9^2)).*(f.^2+12194^2)))+2;
64
+
% Dbruto = Dbruto + ponderacionA;
65
+
% Ebruto = Ebruto + ponderacionA;
66
+
% Lbruto = Lbruto + ponderacionA;
67
+
68
+
% Suma de niveles para los dos rangos
69
+
Dneto(i) = 10*log10(sum(10.^(Dbruto(:)/10)));
70
+
Eneto(i) = 10*log10(sum(10.^(Ebruto(:)/10)));
71
+
Lneto(i) = 10*log10(sum(10.^(Lbruto(:)/10)));
72
+
end
73
+
74
+
% Campo util es la suma del campo directo y el early
75
+
Uneto = 10*log10(10.^(Dneto/10)+10.^(Eneto/10));
76
+
77
+
% Temperatura
78
+
temp = handles.SPLm{handles.Index,17};
79
+
% Velocidad del sonido
80
+
c = 331.4 + 0.6*temp;
81
+
% Densidad del aire. Presion Atm / (constante gas*Temp
82
+
% absoluta)
83
+
rho = 101325 / (287.058*(273.15+temp));
84
+
% Impedancia acustica del aire
85
+
Z = rho*c;
86
+
% Nivel de presion total de la fuente a 1 metro
87
+
SPL = 10*log10(sum(10.^(handles.SPLm{handles.Index,13}/10)));
88
+
% Nivel de potencia de la fuente (SPL a 1 metro). Se supone Q=1
89
+
SWL = SPL - 10*log10(1/(4*pi*1));
90
+
% Potencia de la fuente en vatios
91
+
W = 10.^(SWL/10)*10^-12;
92
+
93
+
%% Representacion grafica del campo Util y el Perjudicial
94
+
95
+
%% Puntos y curva de 0ms a Valor
96
+
% Ordenar valores de distancia de menor a mayor manteniendo su valor de SPl
97
+
% asociado
98
+
Dist_U = sortrows([Distancia',Uneto']);
99
+
% Si tiene distancias duplicadas/iguales se promedia el nivel de las dos
100
+
% distancias iguales y elimina el duplicado
101
+
[C,~,idx] = unique(Dist_U(:,1),'stable');
102
+
val = accumarray(idx,Dist_U(:,2),[],@(x) 10*log10(mean(10.^(x/10))));
103
+
Dist_U = [C,val];
104
+
% A�adir por extrapolacion valor en la posicion de 0 metros hasta la posicion del
105
+
% receptor mas cercano
106
+
%V0_0toVal=interp1(Dist_U(:,1),Dist_U(:,2),0,'linear','extrap');
107
+
%Dist_U = [[0;V0_0toVal]';Dist_U];
108
+
% Separar los valores para graficar y convertirlo en vector fila
109
+
Dist = Dist_U(:,1);
110
+
Urec = Dist_U(:,2);
111
+
% Crear curva potencial. Obtiene coeficientes (Ucoef) y datos de ajuste
112
+
% (Ures)
113
+
[Ucoef,Ures,Uoutput]=fit(Dist,Urec,'power1');
114
+
Uplot=Ucoef.a*Distplot.^Ucoef.b;
115
+
116
+
117
+
%% Curva del campo directo
118
+
% Ordenar valores de distancia de menor a mayor manteniendo su valor de SPl
119
+
% asociado
120
+
Dist_D = sortrows([Distancia',Dneto']);
121
+
% Si tiene distancias duplicadas/iguales se promedia el nivel de las dos
122
+
% distancias iguales y elimina el duplicado
123
+
[C,~,idx] = unique(Dist_D(:,1),'stable');
124
+
val = accumarray(idx,Dist_D(:,2),[],@(x) 10*log10(mean(10.^(x/10))));
125
+
Dist_D = [C,val];
126
+
% A�adir por extrapolacion valor en la posicion de 0 metros hasta la posicion del
127
+
% receptor mas cercano
128
+
%V0_0toVal=interp1(Dist_U(:,1),Dist_U(:,2),0,'linear','extrap');
129
+
%Dist_U = [[0;V0_0toVal]';Dist_U];
130
+
% Separar los valores para graficar y convertirlo en vector fila
131
+
Dist = Dist_D(:,1);
132
+
Drec = Dist_D(:,2);
133
+
% Crear curva potencial. Obtiene coeficientes (Ucoef) y datos de ajuste
134
+
% (Ures)
135
+
[Dcoef,Dres]=fit(Dist,Drec,'power1');
136
+
Dplot=Dcoef.a*Distplot.^Dcoef.b;
137
+
138
+
%% Puntos y curva de valor a infinito ms
139
+
% Ordenar valores de distancia de menor a mayor manteniendo su valor de SPl
140
+
% asociado
141
+
Dist_P = sortrows([Distancia',Lneto']);
142
+
% Si tiene distancias duplicadas/iguales se promedia el nivel de las dos
143
+
% distancias iguales y elimina el duplicado
144
+
[C,~,idx] = unique(Dist_P(:,1),'stable');
145
+
val = accumarray(idx,Dist_P(:,2),[],@(x) 10*log10(mean(10.^(x/10))));
146
+
Dist_P = [C,val];
147
+
% A�adir por extrapolacion valor en la posicion de 0 metros hasta la posicion del
148
+
% receptor mas cercano
149
+
%V0_0toVal=interp1(Dist_P(:,1),Dist_P(:,2),0,'linear','extrap');
150
+
%Dist_P = [[0;V0_0toVal]';Dist_P];
151
+
% Separar los valores para graficar y convertirlo en vector fila
152
+
Dist = Dist_P(:,1);
153
+
Prec = Dist_P(:,2);
154
+
% Obtener coeificentes curva polinomica de grado 1
155
+
[Pcoef,Pres,Poutput] = fit(Dist,Prec,'poly1');
156
+
% Crea la curva
157
+
Pplot = Pcoef.p1*Distplot+Pcoef.p2;
158
+
159
+
160
+
% Claridad EASE
161
+
% Muestra los puntos
162
+
plot(Dist,10*log10( (10.^(Urec/10) + 10.^(Prec/10)) ./ 10^((Dcoef.a*10^Dcoef.b)/10) ) ,'^b','MarkerFaceColor','b'),hold on
163
+
% Muestra la curva
164
+
Curv(1)=plot(Distplot,10*log10( (10.^(Uplot/10) + 10.^(Pplot/10)) ./ 10^((Dcoef.a*10^Dcoef.b)/10) ),'b','LineWidth',1.2);
165
+
166
+
% A�ade las rejillas de las graficas
167
+
grid on
168
+
grid minor
169
+
170
+
% Leyenda
171
+
leyenda{1} = 'G - EASE';
172
+
173
+
%% Teoria corregida
174
+
% Si se ha marcado el calculo de la teoria revisada corregida se calcula y
175
+
% se a�ade a la grafica
176
+
if get(handles.teoriacorregida,'value')==1
177
+
%%% PARAMETROS
178
+
% Nivel de presion total de la fuente a 1 metro
179
+
SPL = 10*log10(sum(10.^(handles.SPLm{handles.Index,13}/10)));
180
+
% Volumen (m3)
181
+
V = handles.SPLm{handles.Index,15};
182
+
% Superficie (m2)
183
+
S = handles.SPLm{handles.Index,16};
184
+
% Mean Free Path
185
+
mfp = 4*V/S;
186
+
% Absorcion media
187
+
alpha = mean(handles.SPLm{handles.Index,11});
188
+
% Absorci�n equivalente de Eyring
189
+
A = -S*log(1-alpha);
190
+
% Temperatura de la sala
191
+
temp = handles.SPLm{handles.Index,17};
192
+
% Velocidad del sonido
193
+
c = 331.4 + 0.6*temp;
194
+
% Densidad del aire. Presion Atm / (constante gas*Temp
195
+
% absoluta)
196
+
rho = 101325 / (287.058*(273.15+temp));
197
+
% Impedancia acustica del aire
198
+
Z = rho*c;
199
+
% Presion de referencia
200
+
pref = 2*10^(-5);
201
+
% Factor de conversion
202
+
Fac = Z/(pref^2);
203
+
% Directividad de la fuente
204
+
Q = str2double(get(handles.Qdir,'string'));
205
+
% Nivel de potencia de la fuente (SPL a 1 metro). Se supone Q=1
206
+
SWL = SPL - 10*log10(1/(4*pi*1));
207
+
% Potencia de la fuente en vatios
208
+
W = 10.^(SWL/10)*10^-12;
209
+
% Tiempo de reverberacion de Eyring
210
+
T = 6*log(10)/c * mfp * 1/(-log(1-alpha));
211
+
% Rango de integracion temporal (s)
212
+
t = handles.Rango * 10^(-3);
213
+
214
+
%% Obtencion de la curvas directo y early
215
+
%% Curva del campo directo
216
+
% Ordenar valores de distancia de menor a mayor manteniendo su valor de SPl
217
+
% asociado
218
+
Dist_D = sortrows([Distancia',Dneto']);
219
+
% Si tiene distancias duplicadas/iguales se promedia el nivel de las dos
220
+
% distancias iguales y elimina el duplicado
221
+
[C,~,idx] = unique(Dist_D(:,1),'stable');
222
+
val = accumarray(idx,Dist_D(:,2),[],@(x) 10*log10(mean(10.^(x/10))));
223
+
Dist_D = [C,val];
224
+
% A�adir por extrapolacion valor en la posicion de 0 metros hasta la posicion del
225
+
% receptor mas cercano
226
+
%V0_0toVal=interp1(Dist_U(:,1),Dist_U(:,2),0,'linear','extrap');
227
+
%Dist_U = [[0;V0_0toVal]';Dist_U];
228
+
% Separar los valores para graficar y convertirlo en vector fila
229
+
Dist = Dist_D(:,1);
230
+
Drec = Dist_D(:,2);
231
+
% Crear curva potencial. Obtiene coeficientes (Ucoef) y datos de ajuste
232
+
% (Ures)
233
+
[Dcoef,Dres]=fit(Dist,Drec,'power1');
234
+
Dplot=Dcoef.a*Distplot.^Dcoef.b;
235
+
236
+
%% Curva del campo early
237
+
% Ordenar valores de distancia de menor a mayor manteniendo su valor de SPl
238
+
% asociado
239
+
Dist_E = sortrows([Distancia',Eneto']);
240
+
% Si tiene distancias duplicadas/iguales se promedia el nivel de las dos
241
+
% distancias iguales y elimina el duplicado
242
+
[C,~,idx] = unique(Dist_E(:,1),'stable');
243
+
val = accumarray(idx,Dist_E(:,2),[],@(x) 10*log10(mean(10.^(x/10))));
244
+
Dist_E = [C,val];
245
+
% A�adir por extrapolacion valor en la posicion de 0 metros hasta la posicion del
246
+
% receptor mas cercano
247
+
%V0_0toVal=interp1(Dist_U(:,1),Dist_U(:,2),0,'linear','extrap');
248
+
%Dist_U = [[0;V0_0toVal]';Dist_U];
249
+
% Separar los valores para graficar y convertirlo en vector fila
250
+
Dist = Dist_E(:,1);
251
+
Erec = Dist_E(:,2);
252
+
% Crear curva potencial. Obtiene coeficientes (Ucoef) y datos de ajuste
253
+
% (Ures)
254
+
[Ecoef,Eres]=fit(Dist,Erec,'power1');
255
+
Eplot=Ecoef.a*Distplot.^Ecoef.b;
256
+
257
+
% Para no mostrar errores cuando se repite el calculo de coeficientes
258
+
warning ('off','all');
259
+
260
+
%% Calculo de coeficiente para el campo directo
261
+
for ii = 1:str2double(get(handles.numIntentos,'string')) % 100 intentos de obtener los coeficientes
262
+
try % Intenta el calculo con puntos de inicio random
263
+
opD = fitoptions('Method','NonlinearLeastSquares',...
264
+
'Robust','LAR',...
265
+
'Lower',[0,0],... % Minimo valor de los coeficientes
266
+
'Upper',[20,20],... % Maximo valor de los coeficientes
267
+
'MaxFunEvals',3000,...
268
+
'MaxIter',3000,...
269
+
'TolFun',10^-3);
270
+
fD = fittype('Fac * W*Q./(4*pi*dist.^2)*Cd',...
271
+
'dependent',{'y'},...
272
+
'independent',{'dist'},...
273
+
'problem', {'Q','W','Fac'},...
274
+
'coefficients',{'Cd'},...
275
+
'options',opD);
276
+
[fitiD,gofD,outputD] = fit(Distplot',10.^(Dplot'/10),fD,'problem',{Q,W,Fac});
277
+
break;
278
+
catch % Si no se ha podido obtener unos coeficientes ejecuta lo siguiente
279
+
if ii==str2double(get(handles.numIntentos,'string'))
280
+
hold off
281
+
ed = errordlg('No se ha podido obtener coeficientes (Directo), intentalo de nuevo','Error');
282
+
set(ed, 'WindowStyle', 'modal');
283
+
uiwait(ed);
284
+
return;
285
+
end
286
+
end
287
+
end
288
+
289
+
290
+
%% Calculo de coeficientes epsilon y Cl para campo late
291
+
for ii = 1:str2double(get(handles.numIntentos,'string')) % 100 intentos de obtener los coeficientes
292
+
try % Intenta el calculo con puntos de inicio random
293
+
fL = fittype('Fac * (4*W)/(S*(-log(1-alpha))) * exp(-(13.82*(dist/c+t)*el/T)) * Cl',...
294
+
'dependent',{'y'},...
295
+
'independent',{'dist'},...
296
+
'problem', {'W','S','alpha','t','T','c','Fac'},...
297
+
'coefficients',{'el','Cl'});
298
+
[fitiL,gofL,outputL] = fit(Distplot',10.^(Pplot'/10),fL,'problem',{W,S,alpha,t,T,c,Fac});
299
+
break;
300
+
catch % Si no se ha podido obtener unos coeficientes ejecuta lo siguiente
301
+
if ii==str2double(get(handles.numIntentos,'string'))
302
+
hold off
303
+
ed = errordlg('No se ha podido obtener coeficientes (Late), intentalo de nuevo','Error');
304
+
set(ed, 'WindowStyle', 'modal');
305
+
uiwait(ed);
306
+
return;
307
+
end
308
+
end
309
+
end
310
+
311
+
%% Calculo de coeficientes epsilon y Ce para campo early
312
+
for ii = 1:str2double(get(handles.numIntentos,'string')) % 100 intentos de obtener los coeficientes
313
+
try % Intenta el calculo con puntos de inicio random
314
+
opE = fitoptions('Method','NonlinearLeastSquares',...
315
+
'Robust','LAR',...
316
+
'MaxFunEvals',3000,...
317
+
'MaxIter',3000,...
318
+
'TolFun',10^-3);
319
+
fE = fittype('Fac * ((4*W)./(S*(-log(1-alpha))*dist) .* (exp(-(13.82*(dist/c)*ee/T))*Ce - exp(-(13.82*(dist/c+t)*el/T)) * Cl))',...
320
+
'dependent',{'y'},...
321
+
'independent',{'dist'},...
322
+
'problem', {'W','S','alpha','T','t','c','Fac','el','Cl'},...
323
+
'coefficients',{'ee','Ce'},...
324
+
'options',opE);
325
+
[fitiE,gofE,outputE] = fit(Distplot',10.^(Eplot'/10),fE,'problem',{W,S,alpha,T,t,c,Fac,fitiL.el,fitiL.Cl});
326
+
break;
327
+
catch % Si no se ha podido obtener unos coeficientes ejecuta lo siguiente
328
+
if ii==str2double(get(handles.numIntentos,'string'))
329
+
hold off
330
+
ed = errordlg('No se ha podido obtener coeficientes (Early), intentalo de nuevo','Error');
331
+
set(ed, 'WindowStyle', 'modal');
332
+
uiwait(ed);
333
+
return;
334
+
end
335
+
end
336
+
end
337
+
warning ('on','all');
338
+
339
+
PDplot = 10*log10(Fac * ( W*Q./(4*pi*Distplot.^2)*fitiD.Cd));
340
+
Peplot = 10*log10(Fac *(...
341
+
(4*W)./(S*(-log(1-alpha))*Distplot) .* ...
342
+
(exp(-(13.82*(Distplot/c)*fitiE.ee/T))*fitiE.Ce - exp(-(13.82*(Distplot/c+t)*fitiL.el/T))*fitiL.Cl)));
343
+
Pdeplot = 10*log10(10.^(PDplot/10)+10.^(Peplot/10));
344
+
Plateplot = 10*log10(Fac .* (4*W)./(S*(-log(1-alpha))) .* exp(-(13.82.*(Distplot./c+t).*fitiL.el/T)) .* fitiL.Cl);
345
+
346
+
Curv(2) = plot(Distplot,10*log10( (10.^(Pdeplot/10) + 10.^(Plateplot/10)) ./ (Fac * ( W./(4*pi*10.^2)*fitiD.Cd)) ),'-.r','LineWidth',1);
347
+
348
+
% Leyendas
349
+
leyenda{2} = 'G - Te�rico';
350
+
351
+
end
352
+
353
+
lgdw=legend(Curv,leyenda);
354
+
lgdw.FontSize = 11;
355
+
%lgdw.Location = 'northeastoutside';
356
+
% Si se utiliza Matlab 2014b o anterior no a�ade el titulo a la leyenda ya
357
+
% que no es compatible
358
+
if ~verLessThan('matlab','8.7')
359
+
title(lgdw,'Curvas')
360
+
end
361
+
ylabel('G (dB)')
362
+
xlabel('Distancia (m)')
363
+
title({'G';...
364
+
['\color{blue} ','Fuente: ',sprintf(': %s',Fuente(1))]});
365
+
366
+
hold off
367
+
+436
include/representarSPLvsDist2.m
+436
include/representarSPLvsDist2.m
···
1
+
function representarSPLvsDist2(hObject,handles)
2
+
3
+
axes(handles.plotfig)
4
+
Rango = handles.Rango;
5
+
SPLmTot = handles.SPLm(:,19); % Matriz con los rayos sin procesar
6
+
SPLm = handles.SPLm(:,20); % Matriz con valores cada ms y por octavas
7
+
Distancia = cell2mat(handles.SPLm(:,6))'; % Vector de distancias fuente-receptor
8
+
Distplot = min(Distancia):0.01:max(Distancia); % Vector de distancias para las curvas
9
+
Fuente = string(handles.SPLm(:,2)');
10
+
Receptor = string(handles.SPLm(:,3)');
11
+
12
+
13
+
% Comprobar que bandas de frecuencias se han elegido
14
+
Bandas = [];
15
+
fsel = [];
16
+
if handles.frec125.Value==1
17
+
Bandas = [Bandas,1];
18
+
fsel = [fsel,125];
19
+
end
20
+
if handles.frec250.Value==1
21
+
Bandas = [Bandas,2];
22
+
fsel = [fsel,250];
23
+
end
24
+
if handles.frec500.Value==1
25
+
Bandas = [Bandas,3];
26
+
fsel = [fsel,500];
27
+
end
28
+
if handles.frec1000.Value==1
29
+
Bandas = [Bandas,4];
30
+
fsel = [fsel,1000];
31
+
end
32
+
if handles.frec2000.Value==1
33
+
Bandas = [Bandas,5];
34
+
fsel = [fsel,2000];
35
+
end
36
+
if handles.frec4000.Value==1
37
+
Bandas = [Bandas,6];
38
+
fsel = [fsel,4000];
39
+
end
40
+
if handles.frec8000.Value==1
41
+
Bandas = [Bandas,7];
42
+
fsel = [fsel,8000];
43
+
end
44
+
if isempty(Bandas); Bandas = 1:7; end
45
+
if isempty(fsel); fsel = 10^3*(2.^((-3:3))); end
46
+
47
+
%% Procesado de datos
48
+
% Separacion e integracion de los niveles antes y despues del rango
49
+
Dneto = zeros(1,numel(SPLm));
50
+
Eneto = zeros(1,numel(SPLm));
51
+
Lneto = zeros(1,numel(SPLm));
52
+
for i=1:numel(SPLm)
53
+
% Direct - Early - Late EXTRACTION
54
+
Dbruto = SPLm{i}(1,Bandas); % Niveles de 0 a 1 ms
55
+
Ebruto = SPLm{i}(2:Rango,Bandas); % Niveles de 1 ms a tiempo de integracion
56
+
Lbruto = SPLm{i}(Rango+1:end,Bandas); % Niveles de tiempo de integracion a infinito
57
+
58
+
% Ponderacion A
59
+
% f = fsel; % Frecuencias seleccionadas
60
+
% ponderacionA = 20*log10(...
61
+
% 12194^2*f.^4 ./...
62
+
% ((f.^2+20.6^2).*sqrt((f.^2+107.7^2).*(f.^2+737.9^2)).*(f.^2+12194^2)))+2;
63
+
% Dbruto = Dbruto + ponderacionA;
64
+
% Ebruto = Ebruto + ponderacionA;
65
+
% Lbruto = Lbruto + ponderacionA;
66
+
67
+
% Suma de niveles para los dos rangos
68
+
Dneto(i) = 10*log10(sum(10.^(Dbruto(:)/10)));
69
+
Eneto(i) = 10*log10(sum(10.^(Ebruto(:)/10)));
70
+
Lneto(i) = 10*log10(sum(10.^(Lbruto(:)/10)));
71
+
end
72
+
73
+
% Campo util es la suma del campo directo y el early
74
+
Uneto = 10*log10(10.^(Dneto/10)+10.^(Eneto/10));
75
+
76
+
%% Representacion grafica del campo Util y el Perjudicial
77
+
colores = colormap(jet(21));
78
+
% Indices de color
79
+
iCU = 1;
80
+
iCUt = 3;
81
+
iCL = 18;
82
+
iCLt = 21;
83
+
%% Puntos y curva de 0ms a Valor
84
+
% Ordenar valores de distancia de menor a mayor manteniendo su valor de SPl
85
+
% asociado
86
+
Dist_U = sortrows([Distancia',Uneto']);
87
+
% Si tiene distancias duplicadas/iguales se promedia el nivel de las dos
88
+
% distancias iguales y elimina el duplicado
89
+
[C,~,idx] = unique(Dist_U(:,1),'stable');
90
+
val = accumarray(idx,Dist_U(:,2),[],@(x) 10*log10(mean(10.^(x/10))));
91
+
Dist_U = [C,val];
92
+
% A�adir por extrapolacion valor en la posicion de 0 metros hasta la posicion del
93
+
% receptor mas cercano
94
+
%V0_0toVal=interp1(Dist_U(:,1),Dist_U(:,2),0,'linear','extrap');
95
+
%Dist_U = [[0;V0_0toVal]';Dist_U];
96
+
% Separar los valores para graficar y convertirlo en vector fila
97
+
Dist = Dist_U(:,1);
98
+
Urec = Dist_U(:,2);
99
+
% Crear curva potencial. Obtiene coeficientes (Ucoef) y datos de ajuste
100
+
% (Ures)
101
+
[Ucoef,Ures,Uoutput]=fit(Dist,Urec,'power1');
102
+
Uplot=Ucoef.a*Distplot.^Ucoef.b;
103
+
% Muestra los puntos
104
+
plot(Dist,Urec,'s','Color',colores(iCU,:),'MarkerFaceColor',colores(iCU,:))
105
+
hold on
106
+
% Muestra los errores
107
+
Uerr = confint(Ucoef);
108
+
Ucurvarec = Ucoef.a*Dist.^Ucoef.b;
109
+
errorbar(Dist,Ucoef.a*Dist.^Ucoef.b,... % Curva
110
+
(Uerr(2)*Dist.^Uerr(4)) - Ucurvarec,... % Error positivo
111
+
Ucurvarec - (Uerr(1)*Dist.^Uerr(3)),... % Error negativo
112
+
'b','LineStyle','none')
113
+
% Muestra la curva
114
+
Curv(1)=plot(Distplot,Uplot,'Color',colores(iCU,:),'LineWidth',1.2);
115
+
116
+
%% Puntos y curva de valor a infinito ms
117
+
% Ordenar valores de distancia de menor a mayor manteniendo su valor de SPl
118
+
% asociado
119
+
Dist_P = sortrows([Distancia',Lneto']);
120
+
% Si tiene distancias duplicadas/iguales se promedia el nivel de las dos
121
+
% distancias iguales y elimina el duplicado
122
+
[C,~,idx] = unique(Dist_P(:,1),'stable');
123
+
val = accumarray(idx,Dist_P(:,2),[],@(x) 10*log10(mean(10.^(x/10))));
124
+
Dist_P = [C,val];
125
+
% A�adir por extrapolacion valor en la posicion de 0 metros hasta la posicion del
126
+
% receptor mas cercano
127
+
%V0_0toVal=interp1(Dist_P(:,1),Dist_P(:,2),0,'linear','extrap');
128
+
%Dist_P = [[0;V0_0toVal]';Dist_P];
129
+
% Separar los valores para graficar y convertirlo en vector fila
130
+
Dist = Dist_P(:,1);
131
+
Prec = Dist_P(:,2);
132
+
% Obtener coeificentes curva polinomica de grado 1
133
+
[Pcoef,Pres,Poutput] = fit(Dist,Prec,'poly1');
134
+
% Crea la curva
135
+
Pplot = Pcoef.p1*Distplot+Pcoef.p2;
136
+
% Muestra los puntos
137
+
plot(Dist,Prec,'^','Color',colores(iCL,:),'MarkerFaceColor',colores(iCL,:))
138
+
% Muestra los errores
139
+
Perr = confint(Pcoef);
140
+
Pcurvarec = Pcoef.p1*Dist+Pcoef.p2;
141
+
errorbar(Dist,Pcoef.p1*Dist+Pcoef.p2,... % Curva
142
+
(Perr(2)*Dist+Perr(4)) - Pcurvarec,... % Error positivo
143
+
Pcurvarec - (Perr(1)*Dist+Perr(3)),... % Error negativo
144
+
'r','LineStyle','none')
145
+
% Muestra la curva
146
+
Curv(2)=plot(Distplot,Pplot,'Color',colores(iCL,:),'LineWidth',1.2);
147
+
148
+
% A�ade las rejillas de las graficas
149
+
grid on
150
+
grid minor
151
+
152
+
% Obtener punto de cruce
153
+
CortesInd = find(abs(Uplot-Pplot)<=(0.01));
154
+
if ~isnan(CortesInd)
155
+
DistCorte = Distplot(CortesInd(1));
156
+
text(DistCorte,Uplot(CortesInd(1))+1,...
157
+
{'EASE',strcat('\bf\color{black}',[sprintf('%4.3f',DistCorte),' m'])},...
158
+
'VerticalAlignment','bottom','HorizontalAlignment','left')
159
+
end
160
+
161
+
%% Guarda un .dat para generar graficas en el LaTex
162
+
%guardarErrores(Ucoef,Pcoef,Distancia,Dist_P,Dist_U)
163
+
%%
164
+
165
+
% Obtener numero de ceros decimales de las pendientes, sumarle 1 y convertir a
166
+
% string, para mostrar tantos numeros decimales en las funciones que se ven en
167
+
% la leyenda
168
+
P11 = num2str(fix(abs(log10(abs(Ucoef.b))))+2);
169
+
P21 = num2str(fix(abs(log10(abs(Pcoef.p1))))+2);
170
+
171
+
% Colores en formato texto para el color de las funciones
172
+
ColU = [num2str(colores(iCU,1)),',',num2str(colores(iCU,2)),',',num2str(colores(iCU,3))];
173
+
ColL = [num2str(colores(iCL,1)),',',num2str(colores(iCL,2)),',',num2str(colores(iCL,3))];
174
+
175
+
% Leyenda
176
+
leyenda{1} = sprintf(['\\bf0 %s %d ms\\rm R^2_{adj} = %4.2f\n\\color[rgb]{%s}y = %4.2f�x^{%4.',P11,'f} \n \\color{white}.'],'a',Rango,Ures.adjrsquare,ColU,Ucoef.a,Ucoef.b);
177
+
leyenda{2} = sprintf(['\\bf%d ms %s\\rm R^2_{adj} = %4.2f\n\\color[rgb]{%s}y = %4.',P21,'f�x+%4.2f \n \\color{white}.'],Rango,'a infinito',Pres.adjrsquare,ColL,Pcoef.p1,Pcoef.p2);
178
+
179
+
%% Teoria corregida
180
+
% Si se ha marcado el calculo de la teoria revisada corregida se calcula y
181
+
% se a�ade a la grafica
182
+
if get(handles.teoriacorregida,'value')==1
183
+
%%% PARAMETROS
184
+
% Nivel de presion total de la fuente a 1 metro
185
+
SPL = 10*log10(sum(10.^(handles.SPLm{handles.Index,13}/10)));
186
+
% Volumen (m3)
187
+
V = handles.SPLm{handles.Index,15};
188
+
% Superficie (m2)
189
+
S = handles.SPLm{handles.Index,16};
190
+
% Mean Free Path
191
+
mfp = 4*V/S;
192
+
% Absorcion media
193
+
alpha = mean(handles.SPLm{handles.Index,11});
194
+
% Absorci�n equivalente de Eyring
195
+
A = -S*log(1-alpha);
196
+
% Temperatura de la sala
197
+
temp = handles.SPLm{handles.Index,17};
198
+
% Velocidad del sonido
199
+
c = 331.4 + 0.6*temp;
200
+
% Densidad del aire. Presion Atm / (constante gas*Temp
201
+
% absoluta)
202
+
rho = 101325 / (287.058*(273.15+temp));
203
+
% Impedancia acustica del aire
204
+
Z = rho*c;
205
+
% Presion de referencia
206
+
pref = 2*10^(-5);
207
+
% Factor de conversion
208
+
Fac = Z/(pref^2);
209
+
% Directividad de la fuente
210
+
Q = str2double(get(handles.Qdir,'string'));
211
+
% Nivel de potencia de la fuente (SPL a 1 metro). Se supone Q=1
212
+
SWL = SPL - 10*log10(1/(4*pi*1));
213
+
% Potencia de la fuente en vatios
214
+
W = 10.^(SWL/10)*10^-12;
215
+
% Tiempo de reverberacion de Eyring
216
+
T = 6*log(10)/c * mfp * 1/(-log(1-alpha));
217
+
% Rango de integracion temporal (s)
218
+
t = handles.Rango * 10^(-3);
219
+
220
+
%% Obtencion de la curvas directo y early
221
+
%% Curva del campo directo
222
+
% Ordenar valores de distancia de menor a mayor manteniendo su valor de SPl
223
+
% asociado
224
+
Dist_D = sortrows([Distancia',Dneto']);
225
+
% Si tiene distancias duplicadas/iguales se promedia el nivel de las dos
226
+
% distancias iguales y elimina el duplicado
227
+
[C,~,idx] = unique(Dist_D(:,1),'stable');
228
+
val = accumarray(idx,Dist_D(:,2),[],@(x) 10*log10(mean(10.^(x/10))));
229
+
Dist_D = [C,val];
230
+
% A�adir por extrapolacion valor en la posicion de 0 metros hasta la posicion del
231
+
% receptor mas cercano
232
+
%V0_0toVal=interp1(Dist_U(:,1),Dist_U(:,2),0,'linear','extrap');
233
+
%Dist_U = [[0;V0_0toVal]';Dist_U];
234
+
% Separar los valores para graficar y convertirlo en vector fila
235
+
Dist = Dist_D(:,1);
236
+
Drec = Dist_D(:,2);
237
+
% Crear curva potencial. Obtiene coeficientes (Ucoef) y datos de ajuste
238
+
% (Ures)
239
+
[Dcoef,Dres]=fit(Dist,Drec,'power1');
240
+
Dplot=Dcoef.a*Distplot.^Dcoef.b;
241
+
242
+
%% Curva del campo early
243
+
% Ordenar valores de distancia de menor a mayor manteniendo su valor de SPl
244
+
% asociado
245
+
Dist_E = sortrows([Distancia',Eneto']);
246
+
% Si tiene distancias duplicadas/iguales se promedia el nivel de las dos
247
+
% distancias iguales y elimina el duplicado
248
+
[C,~,idx] = unique(Dist_E(:,1),'stable');
249
+
val = accumarray(idx,Dist_E(:,2),[],@(x) 10*log10(mean(10.^(x/10))));
250
+
Dist_E = [C,val];
251
+
% A�adir por extrapolacion valor en la posicion de 0 metros hasta la posicion del
252
+
% receptor mas cercano
253
+
%V0_0toVal=interp1(Dist_U(:,1),Dist_U(:,2),0,'linear','extrap');
254
+
%Dist_U = [[0;V0_0toVal]';Dist_U];
255
+
% Separar los valores para graficar y convertirlo en vector fila
256
+
Dist = Dist_E(:,1);
257
+
Erec = Dist_E(:,2);
258
+
% Crear curva potencial. Obtiene coeficientes (Ucoef) y datos de ajuste
259
+
% (Ures)
260
+
[Ecoef,Eres]=fit(Dist,Erec,'power1');
261
+
Eplot=Ecoef.a*Distplot.^Ecoef.b;
262
+
263
+
% Para no mostrar errores cuando se repite el calculo de coeficientes
264
+
warning ('off','all');
265
+
266
+
%% Calculo de coeficiente para el campo directo
267
+
for ii = 1:str2double(get(handles.numIntentos,'string')) % 100 intentos de obtener los coeficientes
268
+
try % Intenta el calculo con puntos de inicio random
269
+
opD = fitoptions('Method','NonlinearLeastSquares',...
270
+
'Robust','LAR',...
271
+
'Lower',[0,0],... % Minimo valor de los coeficientes
272
+
'Upper',[20,20],... % Maximo valor de los coeficientes
273
+
'MaxFunEvals',3000,...
274
+
'MaxIter',3000,...
275
+
'TolFun',10^-3);
276
+
fD = fittype('Fac * W*Q./(4*pi*dist.^2)*Cd',...
277
+
'dependent',{'y'},...
278
+
'independent',{'dist'},...
279
+
'problem', {'Q','W','Fac'},...
280
+
'coefficients',{'Cd'},...
281
+
'options',opD);
282
+
[fitiD,gofD,outputD] = fit(Distplot',10.^(Dplot'/10),fD,'problem',{Q,W,Fac});
283
+
break;
284
+
catch % Si no se ha podido obtener unos coeficientes ejecuta lo siguiente
285
+
if ii==str2double(get(handles.numIntentos,'string'))
286
+
hold off
287
+
ed = errordlg('No se ha podido obtener coeficientes (Directo), intentalo de nuevo','Error');
288
+
set(ed, 'WindowStyle', 'modal');
289
+
uiwait(ed);
290
+
return;
291
+
end
292
+
end
293
+
end
294
+
295
+
296
+
%% Calculo de coeficientes epsilon y Cl para campo late
297
+
for ii = 1:str2double(get(handles.numIntentos,'string')) % 100 intentos de obtener los coeficientes
298
+
try % Intenta el calculo con puntos de inicio random
299
+
fL = fittype('Fac * (4*W)/(S*(-log(1-alpha))) * exp(-(13.82*(dist/c+t)*el/T)) * Cl',...
300
+
'dependent',{'y'},...
301
+
'independent',{'dist'},...
302
+
'problem', {'W','S','alpha','t','T','c','Fac'},...
303
+
'coefficients',{'el','Cl'});
304
+
[fitiL,gofL,outputL] = fit(Distplot',10.^(Pplot'/10),fL,'problem',{W,S,alpha,t,T,c,Fac});
305
+
break;
306
+
catch % Si no se ha podido obtener unos coeficientes ejecuta lo siguiente
307
+
if ii==str2double(get(handles.numIntentos,'string'))
308
+
hold off
309
+
ed = errordlg('No se ha podido obtener coeficientes (Late), intentalo de nuevo','Error');
310
+
set(ed, 'WindowStyle', 'modal');
311
+
uiwait(ed);
312
+
return;
313
+
end
314
+
end
315
+
end
316
+
317
+
%% Calculo de coeficientes epsilon y Ce para campo early
318
+
for ii = 1:str2double(get(handles.numIntentos,'string')) % 100 intentos de obtener los coeficientes
319
+
try % Intenta el calculo con puntos de inicio random
320
+
opE = fitoptions('Method','NonlinearLeastSquares',...
321
+
'Robust','LAR',...
322
+
'MaxFunEvals',3000,...
323
+
'MaxIter',3000,...
324
+
'TolFun',10^-3);
325
+
fE = fittype('Fac * ((4*W)./(S*(-log(1-alpha))*dist) .* (exp(-(13.82*(dist/c)*ee/T))*Ce - exp(-(13.82*(dist/c+t)*el/T)) * Cl))',...
326
+
'dependent',{'y'},...
327
+
'independent',{'dist'},...
328
+
'problem', {'W','S','alpha','T','t','c','Fac','el','Cl'},...
329
+
'coefficients',{'ee','Ce'},...
330
+
'options',opE);
331
+
[fitiE,gofE,outputE] = fit(Distplot',10.^(Eplot'/10),fE,'problem',{W,S,alpha,T,t,c,Fac,fitiL.el,fitiL.Cl});
332
+
break;
333
+
catch % Si no se ha podido obtener unos coeficientes ejecuta lo siguiente
334
+
if ii==str2double(get(handles.numIntentos,'string'))
335
+
hold off
336
+
ed = errordlg('No se ha podido obtener coeficientes (Early), intentalo de nuevo','Error');
337
+
set(ed, 'WindowStyle', 'modal');
338
+
uiwait(ed);
339
+
return;
340
+
end
341
+
end
342
+
end
343
+
warning ('on','all');
344
+
345
+
PDplot = 10*log10(Fac * ( W*Q./(4*pi*Distplot.^2)*fitiD.Cd));
346
+
Peplot = 10*log10(Fac *(...
347
+
(4*W)./(S*(-log(1-alpha))*Distplot) .* ...
348
+
(exp(-(13.82*(Distplot/c)*fitiE.ee/T))*fitiE.Ce - exp(-(13.82*(Distplot/c+t)*fitiL.el/T))*fitiL.Cl)));
349
+
Pdeplot = 10*log10(10.^(PDplot/10)+10.^(Peplot/10));
350
+
Plateplot = 10*log10(Fac .* (4*W)./(S*(-log(1-alpha))) .* exp(-(13.82.*(Distplot./c+t).*fitiL.el/T)) .* fitiL.Cl);
351
+
352
+
Curv(3) = plot(Distplot,Pdeplot,'-.','Color',colores(iCUt,:),'LineWidth',1);
353
+
Curv(4) = plot(Distplot,Plateplot,'-.','Color',colores(iCLt,:),'LineWidth',1);
354
+
355
+
% Colores en formato texto para el color de las funciones
356
+
ColUteo = [num2str(colores(iCUt,1)),',',num2str(colores(iCUt,2)),',',num2str(colores(iCUt,3))];
357
+
ColLteo = [num2str(colores(iCLt,1)),',',num2str(colores(iCLt,2)),',',num2str(colores(iCLt,3))];
358
+
359
+
% Leyendas
360
+
leyenda{3} = sprintf('\\bf0 %s %d ms\\rm Teoria revisada corregida\n \\color[rgb]{%s}\\epsilon_E = %4.3f | C_E = %4.3f | R^2_{adj} = %4.2f \n C_D = %4.3f | R^2_{adj} = %4.2f \n \\color{white}.','a',Rango,ColUteo,fitiE.ee,fitiE.Ce,gofE.adjrsquare,fitiD.Cd,gofD.adjrsquare);
361
+
leyenda{4} = sprintf('\\bf%d ms %s\\rm Teoria revisada corregida\n \\color[rgb]{%s}\\epsilon_L = %4.3f | C_L = %4.3f | R^2_{adj} = %4.2f',Rango,'a infinito',ColLteo,fitiL.el,fitiL.Cl,gofL.adjrsquare);
362
+
363
+
364
+
% Obtener punto de cruce
365
+
CortesInd = find(abs(Pdeplot-Plateplot)<=(0.01));
366
+
if ~isnan(CortesInd)
367
+
DistCorteTeo = Distplot(CortesInd(1));
368
+
text(DistCorteTeo,Pdeplot(CortesInd(1))-1,...
369
+
{'Te�rica',strcat('\bf\color{black}',[sprintf('%4.3f',DistCorteTeo),' m'])},...
370
+
'VerticalAlignment','bottom','HorizontalAlignment','right')
371
+
end
372
+
373
+
% Almacenar variables de interes para despues exportarlas dentro del
374
+
% .mat
375
+
376
+
% Coeficientes
377
+
handles.Elate = fitiL.el;
378
+
handles.Clate = fitiL.Cl;
379
+
handles.Eearly = fitiE.ee;
380
+
handles.Cearly = fitiE.Ce;
381
+
handles.Cdirect = fitiD.Cd;
382
+
383
+
% Parametros
384
+
handles.SPL = SPL; % Nivel de presion a 1 metro
385
+
handles.W = W; % Potencia
386
+
handles.Q = Q; % Directividad
387
+
handles.t = t; % Rango temporal en segundos
388
+
handles.T = T; % Tiempo de reverberaci�n de Eyring
389
+
handles.c = c; % Velocidad de sonido
390
+
handles.S = S; % Superficie
391
+
handles.V = V; % Volumen
392
+
handles.alpha = alpha; % Coeficiente de absorcion medio
393
+
handles.Z = Z; % Impedancia del aire
394
+
handles.rho = rho; % Densidad del aire
395
+
handles.A = A; % Absorcion equivalente de Eyring
396
+
397
+
% Curvas
398
+
handles.Distancia = Distplot; % Vector de distancia
399
+
handles.EaseU = Uplot; % Curva del campo util (EASE)
400
+
handles.EaseP = Pplot; % Curva del campo perjudicial (EASE)
401
+
handles.TeoU = Pdeplot; % Curva del campo util (Teorico)
402
+
handles.TeoP = Plateplot; % Curva del campo perjudicial (Teorico)
403
+
handles.EaseD = Dplot; % Curva del campo directo (EASE)
404
+
handles.EaseE = Eplot; % Curva del campo early (EASE)
405
+
handles.EaseL = Pplot; % Curva del campo late (EASE)
406
+
handles.TeoD = PDplot; % Curva del campo directo (Teorico)
407
+
handles.TeoE = Peplot; % Curva del campo early (Teorico)
408
+
handles.TeoL = Plateplot; % Curva del campo late (Teorico)
409
+
410
+
if exist('DistCorte','var')
411
+
handles.DistCorte = DistCorte;
412
+
end
413
+
if exist('DistCorteTeo','var')
414
+
handles.DistCorteTeo = DistCorteTeo;
415
+
end
416
+
% Actualizar variables en la aplicacion
417
+
guidata(hObject, handles);
418
+
419
+
end
420
+
421
+
lgdw=legend(Curv,leyenda);
422
+
lgdw.FontSize = 11;
423
+
%lgdw.Location = 'northeastoutside';
424
+
% Si se utiliza Matlab 2014b o anterior no a�ade el titulo a la leyenda ya
425
+
% que no es compatible
426
+
if ~verLessThan('matlab','8.7')
427
+
title(lgdw,'Curvas')
428
+
end
429
+
ylabel('Nivel de presi�n ac�stica (dB)')
430
+
xlabel('Distancia (m)')
431
+
title({['SPL seg�n rango de tiempo',sprintf(' - %d ms',Rango)];...
432
+
['\color{blue} ','Fuente: ',sprintf(': %s',Fuente(1))]});
433
+
434
+
hold off
435
+
436
+
+421
include/representarSPLvsDist3.m
+421
include/representarSPLvsDist3.m
···
1
+
function representarSPLvsDist3(hObject,handles)
2
+
3
+
axes(handles.plotfig)
4
+
Rango = handles.Rango;
5
+
SPLmTot = handles.SPLm(:,19); % Matriz con los rayos sin procesar
6
+
SPLm = handles.SPLm(:,20); % Matriz con valores cada ms y por octavas
7
+
Distancia = cell2mat(handles.SPLm(:,6))'; % Vector de distancias fuente-receptor
8
+
Distplot = min(Distancia):0.01:max(Distancia); % Vector de distancias para las curvas
9
+
Fuente = string(handles.SPLm(:,2)');
10
+
Receptor = string(handles.SPLm(:,3)');
11
+
12
+
13
+
% Comprobar que bandas de frecuencias se han elegido
14
+
Bandas = [];
15
+
fsel = [];
16
+
if handles.frec125.Value==1
17
+
Bandas = [Bandas,1];
18
+
fsel = [fsel,125];
19
+
end
20
+
if handles.frec250.Value==1
21
+
Bandas = [Bandas,2];
22
+
fsel = [fsel,250];
23
+
end
24
+
if handles.frec500.Value==1
25
+
Bandas = [Bandas,3];
26
+
fsel = [fsel,500];
27
+
end
28
+
if handles.frec1000.Value==1
29
+
Bandas = [Bandas,4];
30
+
fsel = [fsel,1000];
31
+
end
32
+
if handles.frec2000.Value==1
33
+
Bandas = [Bandas,5];
34
+
fsel = [fsel,2000];
35
+
end
36
+
if handles.frec4000.Value==1
37
+
Bandas = [Bandas,6];
38
+
fsel = [fsel,4000];
39
+
end
40
+
if handles.frec8000.Value==1
41
+
Bandas = [Bandas,7];
42
+
fsel = [fsel,8000];
43
+
end
44
+
if isempty(Bandas); Bandas = 1:7; end
45
+
if isempty(fsel); fsel = 10^3*(2.^((-3:3))); end
46
+
47
+
%% Procesado de datos
48
+
% Separacion e integracion de los niveles antes y despues del rango
49
+
Dneto = zeros(1,numel(SPLm));
50
+
Eneto = zeros(1,numel(SPLm));
51
+
Lneto = zeros(1,numel(SPLm));
52
+
for i=1:numel(SPLm)
53
+
% Direct - Early - Late EXTRACTION
54
+
Dbruto = SPLm{i}(1,Bandas); % Niveles de 0 a 2 ms
55
+
Ebruto = SPLm{i}(2:Rango,Bandas); % Niveles de 2 ms a tiempo de integracion
56
+
Lbruto = SPLm{i}(Rango+1:end,Bandas); % Niveles de tiempo de integracion a infinito
57
+
58
+
% Ponderacion A
59
+
% f = fsel; % Frecuencias seleccionadas
60
+
% ponderacionA = 20*log10(...
61
+
% 12194^2*f.^4 ./...
62
+
% ((f.^2+20.6^2).*sqrt((f.^2+107.7^2).*(f.^2+737.9^2)).*(f.^2+12194^2)))+2;
63
+
% Dbruto = Dbruto + ponderacionA;
64
+
% Ebruto = Ebruto + ponderacionA;
65
+
% Lbruto = Lbruto + ponderacionA;
66
+
67
+
% Suma de niveles para los dos rangos
68
+
Dneto(i) = 10*log10(sum(10.^(Dbruto(:)/10)));
69
+
Eneto(i) = 10*log10(sum(10.^(Ebruto(:)/10)));
70
+
Lneto(i) = 10*log10(sum(10.^(Lbruto(:)/10)));
71
+
end
72
+
73
+
74
+
75
+
%% Representacion grafica del campo Directo/Early/Late
76
+
colores = colormap(jet(21));
77
+
% Indices de color
78
+
iCD = 1;
79
+
iCDt = 3;
80
+
iCE = 6;
81
+
iCEt = 7;
82
+
iCL = 15;
83
+
iCLt = 16;
84
+
%% Puntos y curva del campo directo
85
+
% Ordenar valores de distancia de menor a mayor manteniendo su valor de SPl
86
+
% asociado
87
+
Dist_D = sortrows([Distancia',Dneto']);
88
+
% Si tiene distancias duplicadas/iguales se promedia el nivel de las dos
89
+
% distancias iguales y elimina el duplicado
90
+
[C,~,idx] = unique(Dist_D(:,1),'stable');
91
+
val = accumarray(idx,Dist_D(:,2),[],@(x) 10*log10(mean(10.^(x/10))));
92
+
Dist_D = [C,val];
93
+
% A�adir por extrapolacion valor en la posicion de 0 metros hasta la posicion del
94
+
% receptor mas cercano
95
+
%V0_0toVal=interp1(Dist_U(:,1),Dist_U(:,2),0,'linear','extrap');
96
+
%Dist_U = [[0;V0_0toVal]';Dist_U];
97
+
% Separar los valores para graficar y convertirlo en vector fila
98
+
Dist = Dist_D(:,1);
99
+
Drec = Dist_D(:,2);
100
+
% Crear curva potencial. Obtiene coeficientes (Ucoef) y datos de ajuste
101
+
% (Ures)
102
+
[Dcoef,Dres,Doutput]=fit(Dist,Drec,'power1');
103
+
Dplot=Dcoef.a*Distplot.^Dcoef.b;
104
+
% Muestra los puntos
105
+
plot(Dist,Drec,'s','Color',colores(iCD,:),'MarkerFaceColor',colores(iCD,:)),hold on
106
+
% Muestra los errores
107
+
Derr = confint(Dcoef);
108
+
Dcurvarec = Dcoef.a*Dist.^Dcoef.b;
109
+
errorbar(Dist,Dcoef.a*Dist.^Dcoef.b,... % Curva
110
+
(Derr(2)*Dist.^Derr(4)) - Dcurvarec,... % Error positivo
111
+
Dcurvarec - (Derr(1)*Dist.^Derr(3)),... % Error negativo
112
+
'Color',colores(iCD,:),'LineStyle','none')
113
+
hold on
114
+
% Muestra la curva
115
+
Curv(1)=plot(Distplot,Dplot,'Color',colores(iCD,:),'LineWidth',1.2);
116
+
117
+
%% Curva del campo early
118
+
% Ordenar valores de distancia de menor a mayor manteniendo su valor de SPl
119
+
% asociado
120
+
Dist_E = sortrows([Distancia',Eneto']);
121
+
% Si tiene distancias duplicadas/iguales se promedia el nivel de las dos
122
+
% distancias iguales y elimina el duplicado
123
+
[C,~,idx] = unique(Dist_E(:,1),'stable');
124
+
val = accumarray(idx,Dist_E(:,2),[],@(x) 10*log10(mean(10.^(x/10))));
125
+
Dist_E = [C,val];
126
+
% A�adir por extrapolacion valor en la posicion de 0 metros hasta la posicion del
127
+
% receptor mas cercano
128
+
%V0_0toVal=interp1(Dist_U(:,1),Dist_U(:,2),0,'linear','extrap');
129
+
%Dist_U = [[0;V0_0toVal]';Dist_U];
130
+
% Separar los valores para graficar y convertirlo en vector fila
131
+
Dist = Dist_E(:,1);
132
+
Erec = Dist_E(:,2);
133
+
% Crear curva potencial. Obtiene coeficientes (Ucoef) y datos de ajuste
134
+
% (Ures)
135
+
[Ecoef,Eres,Eoutput]=fit(Dist,Erec,'power1');
136
+
Eplot=Ecoef.a*Distplot.^Ecoef.b;
137
+
% Muestra los puntos
138
+
plot(Dist,Erec,'d','Color',colores(iCE,:),'MarkerFaceColor',colores(iCE,:))
139
+
% Muestra los errores
140
+
Eerr = confint(Ecoef);
141
+
Ecurvarec = Ecoef.a*Dist.^Ecoef.b;
142
+
errorbar(Dist,Ecoef.a*Dist.^Ecoef.b,... % Curva
143
+
(Eerr(2)*Dist.^Eerr(4)) - Ecurvarec,... % Error positivo
144
+
Ecurvarec - (Eerr(1)*Dist.^Eerr(3)),... % Error negativo
145
+
'Color',colores(iCE,:),'LineStyle','none')
146
+
% Muestra la curva
147
+
Curv(2)=plot(Distplot,Eplot,'Color',colores(iCE,:),'LineWidth',1.2);
148
+
149
+
%% Puntos y curva de campor Late
150
+
% Ordenar valores de distancia de menor a mayor manteniendo su valor de SPl
151
+
% asociado
152
+
Dist_L = sortrows([Distancia',Lneto']);
153
+
% Si tiene distancias duplicadas/iguales se promedia el nivel de las dos
154
+
% distancias iguales y elimina el duplicado
155
+
[C,~,idx] = unique(Dist_L(:,1),'stable');
156
+
val = accumarray(idx,Dist_L(:,2),[],@(x) 10*log10(mean(10.^(x/10))));
157
+
Dist_L = [C,val];
158
+
% A�adir por extrapolacion valor en la posicion de 0 metros hasta la posicion del
159
+
% receptor mas cercano
160
+
%V0_0toVal=interp1(Dist_P(:,1),Dist_P(:,2),0,'linear','extrap');
161
+
%Dist_P = [[0;V0_0toVal]';Dist_P];
162
+
% Separar los valores para graficar y convertirlo en vector fila
163
+
Dist = Dist_L(:,1);
164
+
Lrec = Dist_L(:,2);
165
+
% Obtener coeificentes curva polinomica de grado 1
166
+
[Lcoef,Lres,Poutput] = fit(Dist,Lrec,'poly1');
167
+
% Crea la curva
168
+
Lplot = Lcoef.p1*Distplot+Lcoef.p2;
169
+
% Muestra los puntos
170
+
plot(Dist,Lrec,'^','Color',colores(iCL,:),'MarkerFaceColor',colores(iCL,:))
171
+
% Muestra los errores
172
+
Lerr = confint(Lcoef);
173
+
Lcurvarec = Lcoef.p1*Dist+Lcoef.p2;
174
+
errorbar(Dist,Lcoef.p1*Dist+Lcoef.p2,... % Curva
175
+
(Lerr(2)*Dist+Lerr(4)) - Lcurvarec,... % Error positivo
176
+
Lcurvarec - (Lerr(1)*Dist+Lerr(3)),... % Error negativo
177
+
'Color',colores(iCL,:),'LineStyle','none')
178
+
% Muestra la curva
179
+
Curv(3)=plot(Distplot,Lplot,'Color',colores(iCL,:),'LineWidth',1.2);
180
+
181
+
% A�ade las rejillas de las graficas
182
+
grid on
183
+
grid minor
184
+
185
+
% % Obtener punto de cruce
186
+
% CortesInd = find(abs(Dplot-Lplot)<=(0.01));
187
+
% if ~isnan(CortesInd)
188
+
% DistCorte = Distplot(CortesInd(1));
189
+
% text(DistCorte,Dplot(CortesInd(1))+1,...
190
+
% {'EASE',strcat('\bf\color{black}',[sprintf('%4.3f',DistCorte),' m'])},...
191
+
% 'VerticalAlignment','bottom','HorizontalAlignment','left')
192
+
% end
193
+
194
+
195
+
% Obtener numero de ceros decimales de las pendientes, sumarle 1 y convertir a
196
+
% string, para mostrar tantos numeros decimales en las funciones que se ven en
197
+
% la leyenda
198
+
P11 = num2str(fix(abs(log10(abs(Dcoef.b))))+2);
199
+
P21 = num2str(fix(abs(log10(abs(Ecoef.b))))+2);
200
+
P31 = num2str(fix(abs(log10(abs(Lcoef.p1))))+2);
201
+
202
+
% Colores en formato texto para el color de las funciones
203
+
ColD = [num2str(colores(iCD,1)),',',num2str(colores(iCD,2)),',',num2str(colores(iCD,3))];
204
+
ColE = [num2str(colores(iCE,1)),',',num2str(colores(iCE,2)),',',num2str(colores(iCE,3))];
205
+
ColL = [num2str(colores(iCL,1)),',',num2str(colores(iCL,2)),',',num2str(colores(iCL,3))];
206
+
207
+
% Leyendas
208
+
leyenda{1} = sprintf(['\\bf0 %s %d ms\\rm R^2_{adj} = %4.2f\n\\color[rgb]{%s}y = %4.2f�x^{%4.',P11,'f} \n \\color{white}.'],'a',1,Dres.adjrsquare,ColD,Dcoef.a,Dcoef.b);
209
+
leyenda{2} = sprintf(['\\bf%d ms %d ms\\rm R^2_{adj} = %4.2f\n\\color[rgb]{%s}y = %4.2f�x^{%4.',P21,'f} \n \\color{white}.'],1,Rango,Eres.adjrsquare,ColE,Ecoef.a,Ecoef.b);
210
+
leyenda{3} = sprintf(['\\bf%d ms %s\\rm R^2_{adj} = %4.2f\n\\color[rgb]{%s}y = %4.',P31,'f�x+%4.2f \n \\color{white}.'],Rango,'a infinito',Lres.adjrsquare,ColL,Lcoef.p1,Lcoef.p2);
211
+
212
+
%% Teoria corregida
213
+
% Si se ha marcado el calculo de la teoria revisada corregida se calcula y
214
+
% se a�ade a la grafica
215
+
if get(handles.teoriacorregida,'value')==1
216
+
%%% PARAMETROS
217
+
% Nivel total de la fuente a 1 metro
218
+
SPL = 10*log10(sum(10.^(handles.SPLm{handles.Index,13}/10)));
219
+
% Volumen (m3)
220
+
V = handles.SPLm{handles.Index,15};
221
+
% Superficie (m2)
222
+
S = handles.SPLm{handles.Index,16};
223
+
% Mean Free Path
224
+
mfp = 4*V/S;
225
+
% Absorcion media
226
+
alpha = mean(handles.SPLm{handles.Index,11});
227
+
% Absorci�n equivalente de Eyring
228
+
A = -S*log(1-alpha);
229
+
% Temperatura de la sala
230
+
temp = handles.SPLm{handles.Index,17};
231
+
% Velocidad del sonido
232
+
c = 331.4 + 0.6*temp;
233
+
% Densidad del aire. Presion Atm / (constante gas*Temp
234
+
% absoluta)
235
+
rho = 101325 / (287.058*(273.15+temp));
236
+
% Impedancia acustica del aire
237
+
Z = rho*c;
238
+
% Presion de referencia
239
+
pref = 2*10^(-5);
240
+
% Factor de conversion
241
+
Fac = Z/(pref^2);
242
+
% Directividad de la fuente
243
+
Q = str2double(get(handles.Qdir,'string'));
244
+
% Nivel de potencia de la fuente (SPL a 1 metro). Se supone que Q=1
245
+
SWL = SPL - 10*log10(1/(4*pi*1));
246
+
% Potencia de la fuente en vatios
247
+
W = 10.^(SWL/10)*10^-12;
248
+
% Tiempo de reverberacion de Eyring
249
+
T = 6*log(10)/c * mfp * 1/(-log(1-alpha));
250
+
% Rango de integracion temporal (s)
251
+
t = handles.Rango * 10^(-3);
252
+
253
+
% Para no mostrar errores cuando se repite el calculo de coeficientes
254
+
warning ('off','all');
255
+
256
+
%% Calculo de coeficiente para el campo directo
257
+
for ii = 1:str2double(get(handles.numIntentos,'string')) % 100 intentos de obtener los coeficientes
258
+
try % Intenta el calculo con puntos de inicio random
259
+
opD = fitoptions('Method','NonlinearLeastSquares',...
260
+
'Robust','LAR',...
261
+
'Lower',[0,0],... % Minimo valor de los coeficientes
262
+
'Upper',[20,20],... % Maximo valor de los coeficientes
263
+
'MaxFunEvals',3000,...
264
+
'MaxIter',3000,...
265
+
'TolFun',10^-3);
266
+
fD = fittype('Fac * W*Q./(4*pi*dist.^2)*Cd',...
267
+
'dependent',{'y'},...
268
+
'independent',{'dist'},...
269
+
'problem', {'Q','W','Fac'},...
270
+
'coefficients',{'Cd'},...
271
+
'options',opD);
272
+
[fitiD,gofD,outputD] = fit(Distplot',10.^(Dplot'/10),fD,'problem',{Q,W,Fac});
273
+
break;
274
+
catch % Si no se ha podido obtener unos coeficientes ejecuta lo siguiente
275
+
if ii==str2double(get(handles.numIntentos,'string'))
276
+
hold off
277
+
ed = errordlg('No se ha podido obtener coeficientes (Directo), intentalo de nuevo','Error');
278
+
set(ed, 'WindowStyle', 'modal');
279
+
uiwait(ed);
280
+
return;
281
+
end
282
+
end
283
+
end
284
+
285
+
286
+
%% Calculo de coeficientes epsilon y Cl para campo late
287
+
for ii = 1:str2double(get(handles.numIntentos,'string')) % 100 intentos de obtener los coeficientes
288
+
try % Intenta el calculo con puntos de inicio random
289
+
fL = fittype('Fac * (4*W)/(S*(-log(1-alpha))) * exp(-(13.82*(dist/c+t)*el/T)) * Cl',...
290
+
'dependent',{'y'},...
291
+
'independent',{'dist'},...
292
+
'problem', {'W','S','alpha','t','T','c','Fac'},...
293
+
'coefficients',{'el','Cl'});
294
+
[fitiL,gofL,outputL] = fit(Distplot',10.^(Lplot'/10),fL,'problem',{W,S,alpha,t,T,c,Fac});
295
+
break;
296
+
catch % Si no se ha podido obtener unos coeficientes ejecuta lo siguiente
297
+
if ii==str2double(get(handles.numIntentos,'string'))
298
+
hold off
299
+
ed = errordlg('No se ha podido obtener coeficientes (Late), intentalo de nuevo','Error');
300
+
set(ed, 'WindowStyle', 'modal');
301
+
uiwait(ed);
302
+
return;
303
+
end
304
+
end
305
+
end
306
+
307
+
%% Calculo de coeficientes epsilon y Ce para campo early
308
+
for ii = 1:str2double(get(handles.numIntentos,'string')) % 100 intentos de obtener los coeficientes
309
+
try % Intenta el calculo con puntos de inicio random
310
+
opE = fitoptions('Method','NonlinearLeastSquares',...
311
+
'Robust','LAR',...
312
+
'MaxFunEvals',3000,...
313
+
'MaxIter',3000,...
314
+
'TolFun',10^-3);
315
+
fE = fittype('Fac * ((4*W)./(S*(-log(1-alpha))*dist) .* (exp(-(13.82*(dist/c)*ee/T))*Ce - exp(-(13.82*(dist/c+t)*el/T)) * Cl))',...
316
+
'dependent',{'y'},...
317
+
'independent',{'dist'},...
318
+
'problem', {'W','S','alpha','T','t','c','Fac','el','Cl'},...
319
+
'coefficients',{'ee','Ce'},...
320
+
'options',opE);
321
+
[fitiE,gofE,outputE] = fit(Distplot',10.^(Eplot'/10),fE,'problem',{W,S,alpha,T,t,c,Fac,fitiL.el,fitiL.Cl});
322
+
break;
323
+
catch % Si no se ha podido obtener unos coeficientes ejecuta lo siguiente
324
+
if ii==str2double(get(handles.numIntentos,'string'))
325
+
hold off
326
+
ed = errordlg('No se ha podido obtener coeficientes (Early), intentalo de nuevo','Error');
327
+
set(ed, 'WindowStyle', 'modal');
328
+
uiwait(ed);
329
+
return;
330
+
end
331
+
end
332
+
end
333
+
warning ('on','all');
334
+
335
+
PDplot = 10*log10(Fac * ( W*Q./(4*pi*Distplot.^2)*fitiD.Cd));
336
+
Peplot = 10*log10(Fac *(...
337
+
(4*W)./(S*(-log(1-alpha))*(Distplot)) .* ...
338
+
(exp(-(13.82*(Distplot/c)*fitiE.ee/T))*fitiE.Ce - exp(-(13.82*(Distplot/c+t)*fitiL.el/T))*fitiL.Cl)));
339
+
Plateplot = 10*log10(Fac .* (4*W)./(S*(-log(1-alpha))) .* exp(-(13.82.*(Distplot./c+t).*fitiL.el/T)) .* fitiL.Cl);
340
+
341
+
Curv(4) = plot(Distplot,PDplot,'-.','Color',colores(iCDt,:),'LineWidth',1);
342
+
Curv(5) = plot(Distplot,Peplot,'-.','Color',colores(iCEt,:),'LineWidth',1);
343
+
Curv(6) = plot(Distplot,Plateplot,'-.','Color',colores(iCLt,:),'LineWidth',1);
344
+
345
+
% Colores en formato texto para el color de las funciones
346
+
ColDteo = [num2str(colores(iCDt,1)),',',num2str(colores(iCDt,2)),',',num2str(colores(iCDt,3))];
347
+
ColEteo = [num2str(colores(iCEt,1)),',',num2str(colores(iCEt,2)),',',num2str(colores(iCEt,3))];
348
+
ColLteo = [num2str(colores(iCLt,1)),',',num2str(colores(iCLt,2)),',',num2str(colores(iCLt,3))];
349
+
350
+
% Leyendas
351
+
leyenda{4} = sprintf('\\bf0 ms %s %d ms\\rm Teoria revisada corregida\n \\color[rgb]{%s}C_D = %4.3f | R^2_{adj} = %4.2f \n \\color{white}.','a',1,ColDteo,fitiD.Cd,gofD.adjrsquare);
352
+
leyenda{5} = sprintf('\\bf%d ms %s %d ms\\rm Teoria revisada corregida\n \\color[rgb]{%s}\\epsilon_E = %4.3f | C_E = %4.3f | R^2_{adj} = %4.2f \n \\color{white}.',1,'a',Rango,ColEteo,fitiE.ee,fitiE.Ce,gofE.adjrsquare);
353
+
leyenda{6} = sprintf('\\bf%d ms %s\\rm Teoria revisada corregida\n \\color[rgb]{%s}\\epsilon_L = %4.3f | C_L = %4.3f | R^2_{adj} = %4.2f',Rango,'a infinito',ColLteo,fitiL.el,fitiL.Cl,gofL.adjrsquare);
354
+
355
+
356
+
% % Obtener punto de cruce
357
+
% f1a = 10*log10(Peplot);
358
+
% f2a = 10*log10(Plateplot);
359
+
% CortesInd = find(abs(f1a-f2a)<=(0.01));
360
+
% if ~isnan(CortesInd)
361
+
% DistCorteTeo = Distplot(CortesInd(1));
362
+
% text(DistCorteTeo,f1a(CortesInd(1))-1,...
363
+
% {'Te�rica',strcat('\bf\color{black}',[sprintf('%4.3f',DistCorteTeo),' m'])},...
364
+
% 'VerticalAlignment','bottom','HorizontalAlignment','right')
365
+
% end
366
+
367
+
% Almacenar variables de interes para despues exportarlas dentro del
368
+
% .mat
369
+
370
+
% Coeficientes
371
+
handles.Elate = fitiL.el;
372
+
handles.Clate = fitiL.Cl;
373
+
handles.Eearly = fitiE.ee;
374
+
handles.Cearly = fitiE.Ce;
375
+
handles.Cdirect = fitiD.Cd;
376
+
377
+
% Parametros
378
+
handles.SPL = SPL; % Nivel de presion a 1 metro
379
+
handles.W = W; % Potencia
380
+
handles.Q = Q; % Directividad
381
+
handles.t = t; % Rango temporal en segundos
382
+
handles.T = T; % Tiempo de reverberaci�n de Eyring
383
+
handles.c = c; % Velocidad de sonido
384
+
handles.S = S; % Superficie
385
+
handles.V = V; % Volumen
386
+
handles.alpha = alpha; % Coeficiente de absorcion medio
387
+
handles.Z = Z; % Impedancia del aire
388
+
handles.rho = rho; % Densidad del aire
389
+
handles.A = A; % Absorcion equivalente de Eyring
390
+
391
+
% Curvas
392
+
handles.Distancia = Distplot; % Vector de distancia
393
+
handles.EaseD = Dplot; % Curva del campo directo (EASE)
394
+
handles.EaseE = Eplot; % Curva del campo early (EASE)
395
+
handles.EaseL = Lplot; % Curva del campo late (EASE)
396
+
handles.TeoD = PDplot; % Curva del campo directo (Teorico)
397
+
handles.TeoE = Peplot; % Curva del campo early (Teorico)
398
+
handles.TeoL = Plateplot; % Curva del campo late (Teorico)
399
+
400
+
% Actualizar variables en la aplicacion
401
+
guidata(hObject, handles);
402
+
403
+
end
404
+
405
+
lgdw=legend(Curv,leyenda);
406
+
lgdw.FontSize = 11;
407
+
lgdw.Location = 'bestoutside';
408
+
%lgdw.Location = 'northeastoutside';
409
+
% Si se utiliza Matlab 2014b o anterior no a�ade el titulo a la leyenda ya
410
+
% que no es compatible
411
+
if ~verLessThan('matlab','8.7')
412
+
title(lgdw,'Curvas')
413
+
end
414
+
ylabel('Nivel de presi�n ac�stica (dB)')
415
+
xlabel('Distancia (m)')
416
+
title({['SPL seg�n rango de tiempo',sprintf(' - %d ms',Rango)];...
417
+
['\color{blue} ','Fuente: ',sprintf(': %s',Fuente(1))]});
418
+
419
+
hold off
420
+
421
+