model to generate Data

forward_model/us_mp_model.m

+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
+