ita_diffraction_plot_utd_total_field_stienen_diss_w_ground.m 4.02 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
%% Total sound field of a difraction problem at a rectagular wedge using UTD

% Author: Jonas Stienen


%% Scene

n1  = [ 0  1  0 ]; % main face 1
n2  = [ 1  0  0 ]; % opposite face 2 / facade
n3  = [ 0  1  0 ]; % ground plane / face 3
loc = [  0  0  0 ];
w = itaInfiniteWedge( n1, n2, loc );

r    = [  2  -2  0 ];
r2   = [ -2  -2  0 ]; % image receiver 2
r3   = [  2  -4  0 ]; % image receiver 3
r23  = [ -2  -4  0 ]; % image receiver 23 /  32
s1  = 5 * [ 0  -1  0 ];

a = w.apex_point( s1, r );
ad = w.edge_direction;
c = 344; % speed of sound

f = [ 20, 50, 100, 200, 400, 800, 1600, 3200, 6400, 12800, 24000 ]';
f_legend = { '20 Hz', '50 Hz', '100 Hz', '200 Hz', '400 Hz', '800 Hz', ...
    '1.6 kHz', '3.2 kHz', '6.4 kHz', '12.8 kHz', '24 kHz' };

%f = [ 20, 50, 100, 200, 400 ]';
%f_legend = { '20 Hz', '50 Hz', '100 Hz', '200 Hz', '400 Hz' };

num_angles = 1000;
alpha_d_start = 0;
alpha_d_end = pi;
alpha_d = linspace( alpha_d_start, alpha_d_end, num_angles );

% Set different receiver positions rotated around the aperture
S = norm( s1 ) * [ - cos( alpha_d ); sin( alpha_d ); zeros( 1, numel( alpha_d ) ) ]';

% wavenumber
k = 2 * pi * f ./ c;

% store results in ita audio format
res_data_template = itaResult();
res_data_template.resultType = 'signal';
res_data_template.freqVector = f;


%% Sound field simulation

N = size( S, 1 );
for n = 1 : N
    
    s = S( n, 1:3 );
    
    % Zero order
    
    r_dir = norm( s - r );
    H_i = 1 ./ r_dir .* exp( -1i .* k .* r_dir ); % direct sound component
    
    r_1_dir = norm( s - r2 );
    H_r = 1 ./ r_1_dir .* exp( -1i .* k .* r_1_dir ); % reflected sound component

    H_d = ita_diffraction_utd( w, s, r, f, c ); % diffracted sound component
    
    % Ground reflection image
    
    r_dir_1 = norm( s - r3 );
    H_i_1 = 1 ./ r_dir_1 .* exp( -1i .* k .* r_dir_1 ); % direct sound component
    
    r_1_dir_1 = norm( s - r23 );
    H_r_1 = 1 ./ r_1_dir_1 .* exp( -1i .* k .* r_1_dir_1 ); % reflected sound component

    [ H_d_1 , D_half ] = ita_diffraction_utd( w, s, r3, f, c ); % diffracted sound component
    
    % zero and first order components
    
    H = H_d; % Total sound field ... diffracted component always present
    if ~ita_diffraction_shadow_zone( w, s, r )
       H = H + H_i;
    end
    if ita_diffraction_reflection_zone( w, s, r, false ) % opposite face reflection
        H = H + H_r;
    end
    
    H_1 = H_d_1;
    if ~ita_diffraction_shadow_zone( w, s, r3 )
        H_1 = H_1 + H_i_1;
    end
    if ita_diffraction_reflection_zone( w, s, r3, false ) % opposite face reflection
        H_1 = H_1 + H_r_1;
    end
    
    H_1_half = zeros( numel( H_d_1 ), 1 );
    if ita_diffraction_shadow_zone( w, s, r3 )
        H_1_half = H_d_1; % Shadow zone yes diffraction only
    else
        H_1_half = H_i_1; % Direct sound yes
        if ~ita_diffraction_reflection_zone( w, s, r3, false )
            if ita_diffraction_utd_illumination_sign( w, s, r3, false ) == -1
                sf = ita_diffraction_utd_illumination_fadeout( w, s, r3, false );
                H_1_half = H_1_half + sf .* H_d_1; % closer to shadow zone diffraction yes
            end
        else
            humbug = 1;
        end
    end
    
    att_sum = res_data_template;
109
    att_sum.freqData = ( H_1_half ) ./ H_i; % Insertion loss
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
        
    if n == 1
        att_sum_all = att_sum;
    else
        att_sum_all = ita_merge( att_sum_all, att_sum );
    end
    
end


%% Plot

fg = figure( 'units', 'normalized', 'outerposition', [0 0 1 1], 'visible', 'on' );

h = plot( rad2deg( alpha_d ), flip( att_sum_all.freqData_dB', 2 ) );

title( 'Transfer function nomalized to free-field transmission' )
xlabel( 'alpha_d / degree' )
ylabel( 'dB' )

lgd = legend( flip( f_legend ), 'Location', 'southwest' );
lgd.NumColumns = 3;

xlim( [ -5 185 ] );
ylim( [ -35 10 ] );
grid on;

lw = linspace( 0.5, 2.5, numel( f ));
for i=1:numel( f )
    h( i ).LineWidth = lw( i );
end


%% Export

set( fg, 'papersize', [ 16*1.5 9*1.5 ], 'units', 'centimeters' )
print( 'ita_diffraction_plot_utd_total_field_stienen_diss_with_ground', '-fillpage', '-dpdf' )