%% 15th February 2021, Louise Wilson %Plot all of the boat transects in terms of pressure and particle %acceleration, alongside boat passes (distance of boat from array). load('ArrayData_11Aug2021_GoatIsland_Koura_+-2_5secs.mat'); %load pressure data %correct datetime in 'gps' as in wrong format for i=1:height(gps); date=datetime(gps.LOCALDATE{i},'InputFormat','yyyy/MM/dd'); time=gps.LOCALTIME(i); gps.LOCALDT(i)=date+time; end %% Plot boat range and speed % Calculate horizontal distances between source and receiver array hydroLat=-36.266833; %hydrophone position hydroLon=174.790417; HydroDepth = 0.0085; %depth of hydrophone (km) SourceDepth = 0.0012 ; %depth of source (km) for j = 1:height(gps); dist=deg2km(distance(gps.LAT(j),gps.LONGITUDE(j),hydroLat,hydroLon)); A = dist;%horizontal distances (km) B = HydroDepth - SourceDepth; %calculate diff. b/w source + receiver depths (km) C = sqrt(A^2+B^2); %slope distance (km) gps.slopeDist(j) = C*1000; %metres gps.horzDist(j)=A*1000; %metres end %Add index for seconds (to sync plot with spectrograms) boat_data=gps(1:296,:); %296 source positions idx=[0:5:1475] %5 second interval between source positions boat_data.idx(:)=idx; %(Note that whilst source positions were recorded every five seconds, there %are sometimes gaps of longer than five seconds (e.g. 10, 15) between %time of each source position. This index allows boat source positions %to be aligned with spectrogram data, where all five seconds chunks are %concatenated (irrespective of gap in timing between each source position) %-using time would not work.) % Plot subplot(3,1,1) set(gca,'layer','top','fontname','arial'); plot(boat_data.idx,boat_data.horzDist) %boat range (horizontal) ylabel('Boat range (m)','fontsize',12); yyaxis right plot(boat_data.idx,boat_data.SPEED) %speed of boat (km) ylabel('Boat speed (km/h)','Rotation',270,'VerticalAlignment','bottom',... 'fontsize',12) set(gca,'tickdir','out'); xticks([0:200:1400]) xlim([0 1480]) clearvars -except gps %% ACCELERATION % Calculate PM quantities from equations in Peter Rogers' .docx. fs = 72000; %sample rate used to record all hydrophone data rho = 1026; % density of water (kg/m^-3) h = 0.4; % distance between hydrophones in metres for z=1:height(gps) %Converting micropascals to pascals pos1_Pa=(gps.pos1{z,1}/power(10,6)); %units Pa pos2_Pa=(gps.pos2{z,1}/power(10,6)); pos3_Pa=(gps.pos3{z,1}/power(10,6)); pos4_Pa=(gps.pos4{z,1}/power(10,6)); pos5_Pa=(gps.pos5{z,1}/power(10,6)); pos6_Pa=(gps.pos6{z,1}/power(10,6)); %Pressure differential approximation gps.Ax{z}=(pos2_Pa-pos1_Pa)/(rho*h); gps.Ay{z}=(pos4_Pa-pos3_Pa)/(rho*h); gps.Az{z}=(pos6_Pa-pos5_Pa)/(rho*h); %Magnitude (m/s/s) Axyz=sqrt((gps.Ax{z}.^2)+(gps.Ay{z}.^2)+(gps.Az{z}.^2)); %Magnitude (µm/s/s) gps.Axyz{z}=Axyz*power(10,6); clearvars -except gps fs rho h n end plot_data=gps(1:296,:); magnitude=cell2mat(plot_data.Axyz); nfft = 2048; % adjust fft size overlap = nfft/2; window = nfft; [S,F,T,P] = spectrogram(magnitude,hamming(window),overlap,nfft,fs,'yaxis'); Fcorrected = 2*pi*F; %convert Eulers equation from velocity to acceleration PP=zeros(length(F), length(T)); for i=1:length(Fcorrected); %Fcorr=Fcorrected'; Pcorrected=(Fcorrected(i))*P(i,:); PP(i,:) = Pcorrected; end %freq_res=fs/nfft; %35.1563 Hz %temp_res=overlap/fs; %0.0142 seconds subplot(3,1,3) surf(T,F,10*log10(PP),'edgecolor','none'); colormap(jet); axis tight; view(0,90); grid off set(gca,'tickdir','out','layer','top','fontname',... %change size of axes 'arial','fontsize',10); xlabel('Time', 'fontsize', 12,'FontName','Arial'); ylabel('Frequency (Hz)','fontsize',12, 'FontName', 'Arial'); ylabel(colorbar,'PSD (dB re (1µm s^{-2})^2 / Hz)','Rotation',270,'VerticalAlignment','bottom',... 'fontsize',12); %positions colorbar axis set(gca,'YScale','log'); %set y-axis scale to logarithmic ylim ([100 3000]); caxis([30 90]); xlabel('Seconds'); clearvars -except gps %% PRESSURE: PSD analysis on mean data per source position (1 Hz bins) %Input variables (same values as for acceleration) fs=72000; %sample rate nfft=2048; %fft length window=nfft; %window size overlap=nfft/2; %amount of overlap between segments (50%) %freq_res=fs/nfft; %35.1563 Hz %temp_res=overlap/fs; %0.0142 seconds plot_data=gps(1:296,:); %subset data to end at 10:24:00 idx=[1:1:296]; plot_data.idx(:)=idx; pressure_data=cell2mat(plot_data.posMean); [S,F,T,P] = spectrogram(pressure_data,hamming(window),overlap,nfft,fs,'yaxis'); %Plot subplot(3,1,2) surf(T,F,10*log10(P),'edgecolor','none'); colormap(jet); axis tight; view(0,90); grid off set(gca,'layer','top','fontname','arial','fontsize',10); set(gca,'YScale','log'); %set y-axis scale to logarithmic ylabel('Frequency (Hz)','fontsize',12); ylim([100 3000]); grid off colorbar caxis([30 120]); %sets scale on colour bar ylabel(colorbar,'PSD (dB re 1µPa^2 / Hz)','Rotation',270,'VerticalAlignment','bottom',... 'fontsize',12); %positions colorbar axis set(gca,'tickdir','out'); %% Save f = gcf; exportgraphics(f,'Supplementary Figure S2.pdf','Resolution',600) save('Supplementary Figure S2.fig')