us_mp_model.m 7.03 KB
Newer Older
John, Fabian's avatar
John, Fabian committed
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
classdef us_mp_model
    %US_MP_MODEL This class implements the multipath propagation ultrasonic
    %fresnel model.
    %   Refer publication: Fabian John, Marco Cimdins, Horst Hellbrueck,
    %   Underwater Ultrasonic Multipath Diffraction Model for 
    %   Communication and Sensing Applications, IEEE-Sensors, 2021
    
    properties
        d1;          % LOS-projection distance obstacle to Tx
        d2;          % LOS-projection distance obstacle to Rx
        b;           % obstacle diameter
        f_low;       % lower frequency bound
        f_high;      % upper frequency bound
        T_pulse;     % pulse length
        d_refl;      % distances vector for reflections
        R;           % reflection factor vector
        c = 1503;    % underwater sound velocity
        F0_los = 1;  % LOS transmission factor
        f_res = 100; % frequency resolution factor
    end
    
    methods
        function obj = us_mp_model(d1, d2, b, f_low, f_high, T_pulse, d_refl, R)
            %US_MP_MODEL Construct an instance of this class
            obj.d1 = d1;
            obj.d2 = d2;
            obj.b = b;
            obj.f_low = f_low;
            obj.f_high = f_high;
            obj.T_pulse = T_pulse;
            obj.d_refl = d_refl;
            obj.R = R;
        end
        
        function [f, y_fft] = get_rx_db(obj,d_los)
            %METHOD1 returns the spectrum for a given obstacle distance
            %perpendicular to the los
            [t_front, y_front] = obj.get_diffr_front(d_los);
            [~, y_back] = obj.get_diffr_back(d_los);
            [~, y_lost] = obj.get_los_transmission(d_los);
            [~, y_reflC] = obj.get_refl();
            
            tar_length = 2*length(t_front);
            t_tar = [t_front, t_front+t_front(end)-t_front(1)+t_front(2)-t_front(1)];
            y_front = obj.append_zeros(y_front, tar_length);
            y_back = obj.append_zeros(y_back, tar_length);
            y_lost = obj.append_zeros(y_lost, tar_length);
            y_refl = zeros(length(y_reflC), tar_length);
            for idx = 1:length(y_reflC)
                y_refl(idx,:) = obj.append_zeros(y_reflC{idx}, tar_length);
            end
            
            y_total = y_front + y_back + y_lost + sum(y_refl, 1);
            [y_fft, f] = fft_single_sided(t_tar, y_total, obj.f_high);
            y_fft = 10*log10(y_fft/sqrt(2))-10*log10(1+sum(obj.R));
        end
    end
    
    methods (Access = private)
        function f = get_mod_frequencies(obj)
            f = obj.f_low*.9:(obj.f_high*1.1 - obj.f_low*.9)/obj.f_res:obj.f_high*1.1;
        end
        
        function fres = get_mod_f_resolution(obj)
            f = obj.get_mod_frequencies();
            fres = ((f(2)-f(1))/((1)*length(f)));
        end
        
        function [t, y_t] = get_diffr_front(obj, d_los)
            % calculate the front diffraction component first
            f = obj.get_mod_frequencies();
            lambda = obj.c./f;
            dLosFront = d_los - (obj.b/2);
            v_front = dLosFront.*sqrt( (2*(obj.d1 + obj.d2))./(lambda*obj.d1.*obj.d2) );
            F_front = ((1+1i)/2)*...
            ( (0.5 - fresnelc(-v_front)) - ...
            1i*(0.5*fresnels(-v_front)) );
        
            % now calculate the attenuation in time domain for the front
            % component
            mag = obj.get_mod_f_resolution()*abs(F_front);
            phi = angle(F_front);
            

            t = -15/f(1):obj.get_fs():15/f(1);
            y_t = zeros(1,length(t));
            for idx = 1:length(t)
                y_t(idx) = sum(mag.*sin(2*pi*f*t(idx)+phi), 2);
            end
            y_t = y_t.*obj.get_hanning_zerofilled(t);
        end
        
        function [t, y_t] = get_diffr_back(obj, d_los)
            % calculate the front diffraction component first
            f = obj.get_mod_frequencies();
            lambda = obj.c./f;
            dLosBack = d_los + (obj.b/2);
            v_back = dLosBack.*sqrt( (2*(obj.d1 + obj.d2))./(lambda*obj.d1.*obj.d2) );
            F_back = ((1+1i)/2)*...
            ( (0.5 - fresnelc(v_back)) - ...
            1i*(0.5*fresnels(v_back)) );
            
            
            % now calculate the attenuation in time domain for the front
            % component
            mag = obj.get_mod_f_resolution()*abs(F_back);
            phi = angle(F_back);

            t = -15/f(1):obj.get_fs():15/f(1);
            y_t = zeros(1,length(t));
            for idx = 1:length(t)
                y_t(idx) = sum(mag.*sin(2*pi*f*t(idx)+phi), 2);
            end
            
            delta_d = sqrt(obj.d1^2+dLosBack^2) + sqrt(obj.d2^2+dLosBack^2) - (obj.d1+obj.d2);
            dt_back = delta_d/obj.c;
            alpha = ((obj.d1+obj.d2)/((dt_back*obj.c/2)+obj.d1+obj.d2))^2;
            y_t = y_t.*obj.get_hanning_zerofilled(t)*alpha;
            
            % phase accurate superposition -> shift arrival of back
            % component by additional travel path around the obstacle
            [t, y_t] = obj.phase_correction_shift(t, y_t, dt_back);
        end
        
        function [t, y_t] = get_los_transmission(obj, d_los)
            f = obj.get_mod_frequencies();
            t = -15/f(1):obj.get_fs():15/f(1);
            
            obst_sz = 0.5*obj.b;
            y_t = zeros(1,length(t));
            if abs(d_los) > obst_sz
                return;
            else
                magTm = obj.F0_los*cos(pi/2*abs(d_los)/(obst_sz));
            end
            for idx = 1:length(t)
                y_t(idx) = sum(magTm.*sin(2*pi*f*t(idx)), 2);
            end
            y_t = y_t.*obj.get_hanning_zerofilled(t);
        end
        
        function [tout, yout] = get_refl(obj)
            f = obj.get_mod_frequencies();
            tout = cell(length(obj.R), 1);
            yout = cell(length(obj.R), 1);
            for r_idx = 1:length(obj.R)
                t = -15/f(1):obj.get_fs():15/f(1);
                y_t = zeros(1,length(t));
                for idx = 1:length(t)
                    y_t(idx) = sum(obj.R(r_idx).*sin(2*pi*f*t(idx)), 2);
                end
                y_t = y_t.*obj.get_hanning_zerofilled(t);
                % shift by additional path of reflection distance
                [t, y_t] = obj.phase_correction_shift(t, y_t, obj.d_refl(r_idx)/obj.c);
                tout{r_idx} = t;
                yout{r_idx} = y_t;
            end
        end
        
        function [t, y_t] = phase_correction_shift(obj, t, y_t, dt)
            shift_idx = find(t >= t(1)+dt, 1) - 1;
            y_t = [zeros(1, shift_idx), y_t];
            t = [t, t(end)+obj.get_fs()*(1:shift_idx)];
        end
        
        function y_hann = get_hanning_zerofilled(obj, t)
            h_size = ceil(obj.T_pulse/obj.get_fs());
            y_hann = [zeros(1, ceil((length(t)-h_size)/2)), hanning(h_size)', zeros(1, floor((length(t)-h_size)/2))];
        end
        
        function vec = append_zeros(obj, vec_in, tar_length)
            vec = [vec_in, zeros(1, tar_length - length(vec_in))];
        end
        
        function fs = get_fs(obj)
            f = obj.get_mod_frequencies;
            fs = 1/(f(end)*20);
        end
    end
end