diff --git a/dataProcessing/LTSAs/longTermSpectra_GoMex2018.asv b/dataProcessing/LTSAs/longTermSpectra_GoMex2018.asv new file mode 100644 index 0000000..85c58d6 --- /dev/null +++ b/dataProcessing/LTSAs/longTermSpectra_GoMex2018.asv @@ -0,0 +1,122 @@ +% script modified by J. Haxel May 2015 to calculate the long-term time +% averaged spectra of acoustic signals at NRS +% system response = (preAmp + instResp) is removed and energy values are +% normalized to dB re 1uPa @ 1m +% +% CF2 data - need to use *CF2*.m for all subroutines +% +% data info +%c +if (1) + datadir = 'E:\sg639\wav\sg639-1kHz\'; + locname = 'GoMex'; + fnums = 1:49999; + chan = 1; + frameSize = 1000; + fdates = 'E:\sg639\wav\sg639-1kHz\file_dates-sg639_GoMex_May18-1kHz.txt'; + sRate = 1000; +elseif(0) + +end + % +% Parameters for time averaging and start/stop times +period = 1/86400; %5/1440 %1/24 %5/1440 %hourly =1/24 %600/(60*60*24); % 600 sec in days +baseyear = 'auto'; % auto means get it from first day of data +zeroPad = 0; % zero-pad the signal +nOverlap = 0; % overlap data points +% +%%%%%%%%%%%%%%%%%%%%%%% Configure data Grab %%%%%%%%%%%%%%%%%%%%%%%%%%%% +% +hIndex = readHarufileDateIndex(fdates); +if (0) %change for the old format of datafile + fn0 = sprintf('%sdatafile.%03d', datadir, fnums(1)); + fn1 = sprintf('%sdatafile.%03d', datadir, fnums(end)); +elseif (0) + fn0 = sprintf('%s%08d.DAT', datadir, fnums(1)); + fn1 = sprintf('%s%08d.DAT', datadir, fnums(end)); +else + fn0 = [datadir hIndex.hname{1}(3:end)]; + fn1 = [datadir hIndex.hname{end}(3:end)]; +end +% NOT NEEDED FOR WORKING WITH WAV FILES +% [dtv0,yr,jd,hr,mn,sc,sRate,t,y]=readCF2hdr(fn0,0); +% [dtv1,yr,jd,hr,mn,sc,sRate,t,y]=readCF2hdr(fn1,0); +dtv0 = hIndex.time(1); +dtv1 = hIndex.time(end); +% +if (ischar(baseyear) && strcmp(baseyear, 'auto')) +% baseyear = sub(datevec(dtv0), 1); % new dt encoding from soundIn +baseyear = sub(datevec(hIndex.time(1)), 1); +end +% +% Get start and stop times that include only whole intervals. +dt0 = ceil(dtv0/period) * period; % get start of first whole period +dt1 = floor(dtv1/period) * period; % get end of last whole period +% +% get the system response to subtract out of power vaues later +%[SR] = get_sys_response(PA_rev,gain,filt_sw,sRate,frameSize); +%[SR] = get_sys_response(9,0,4,double(sRate),frameSize); + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% %%%%%%%%%%%%%%%%% Start the loop for calculating LTSpect %%%%%%%%%%%%%% +% +if (1) + gram = zeros((frameSize+zeroPad)/2, round((dt1-dt0)/period)); + fRate = sRate / (frameSize - nOverlap); + printf('%d dots:', nCols(gram)) + for i = 1 : nCols(gram) + fprintf(1,'.%s', iff(mod(i,60) == 0, char(10), '')); + dt = dt0 + (i-1) * period; +% Do around 1 million samples at a time. + t0 = dt; + gSum = zeros((frameSize+zeroPad)/2, 1); + ng = 0; +% while (t0 < dt + period) + %t1 = min([dt+period double(t0+1e6/sRate/24/60/60)]); %for big AVG's + t1 = dt + period; + %xt = haruSoundAtTime_CF2_2channel([t0 t1], hIndex, datadir, chan); +% xt = haruSoundAtTime_CF2([t0 t1], hIndex, datadir); %note - calling CF2 specific routine + xt = haruSoundAtTime([t0 t1], hIndex, datadir); + + %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + %xt_2 = xt{1}; %xt_3 = xt_2(:,1); %extra step do to 2 channel data + %x_mn = xt_3 - mean(xt_3); %have to do this in 2 steps now because 2 channel data + x_mn = xt{1} - mean(xt{1}); %remove the sample mean from the chunk of data for 16bit + x = (2.5/(2^16)).*x_mn; % convert to Volts by (Voltrange/bitrange)- for 16bit + %x = (4/(2^12)).*x_mn; %convert to Volts -12bit + g = joespect_LTacoustic(x, frameSize, nOverlap, zeroPad, 'hann', sRate); %5000 for sample rate + +% [g,f,t] = spectrogram(x,hann(frameSize),nOverlap,frameSize,sRate); +% g = abs(g(2:end)); + % last input on line above is sample rate + if (nCols(g) == 0), break; end % skip partial frame at end of day + gSum = gSum + mean(g,2); % sum up the energy - will divide later + ng = ng + nCols(g); + %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + %t0 = t0 + nCols(g)/.9765624 / 24/60/60; %have to change this to make it run m10 hickup - + t0 = t0 + nCols(g) / fRate / 24/60/60; %change back to this + %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% end + grm = 10*log10(gSum); %to convert power to dB + % + %GRAM = grm - SR'; + % + %gram(:,i) = GRAM; %must add this in to get system response removed + gram(:,i) = grm; + % +% + end + fprintf(1, '\n'); +end +% remove the system gain from the spectra +%LTgram = gram - repmat(SR',1,length(gram(1,:))); +% +%cd E:\ROSS_SEA-2014 +%cd C:\haxel\OCEAN_NOISE\YAQ_HEAD\ +%cd C:\haxel\LAU_BASIN\ANALYSIS\MATAS +%% +cd C:\Users\haver\Documents\NRS\NRS10\LTSA + +save NRS10_20162017_5minAvgSpect_cal070118 gram dt0 dt1 frameSize sRate + +% \ No newline at end of file diff --git a/dataProcessing/LTSAs/longTermSpectra_GoMex2018.m b/dataProcessing/LTSAs/longTermSpectra_GoMex2018.m new file mode 100644 index 0000000..50cc7a5 --- /dev/null +++ b/dataProcessing/LTSAs/longTermSpectra_GoMex2018.m @@ -0,0 +1,117 @@ +% script modified by J. Haxel May 2015 to calculate the long-term time +% averaged spectra of acoustic signals at NRS +% system response = (preAmp + instResp) is removed and energy values are +% normalized to dB re 1uPa @ 1m +% +% CF2 data - need to use *CF2*.m for all subroutines +% +% data info +%c +if (1) + datadir = 'E:\sg639\wav\sg639-1kHz\'; + locname = 'GoMex'; + fnums = 1:49999; + chan = 1; + frameSize = 1000; + fdates = 'E:\sg639\wav\sg639-1kHz\file_dates-sg639_GoMex_May18-1kHz.txt'; + sRate = 1000; +elseif(0) + +end + % +% Parameters for time averaging and start/stop times +period = 1/86400; %5/1440 %1/24 %5/1440 %hourly =1/24 %600/(60*60*24); % 600 sec in days +baseyear = 'auto'; % auto means get it from first day of data +zeroPad = 0; % zero-pad the signal +nOverlap = 0; % overlap data points +% +%%%%%%%%%%%%%%%%%%%%%%% Configure data Grab %%%%%%%%%%%%%%%%%%%%%%%%%%%% +% +hIndex = readHarufileDateIndex(fdates); +if (0) %change for the old format of datafile + fn0 = sprintf('%sdatafile.%03d', datadir, fnums(1)); + fn1 = sprintf('%sdatafile.%03d', datadir, fnums(end)); +elseif (0) + fn0 = sprintf('%s%08d.DAT', datadir, fnums(1)); + fn1 = sprintf('%s%08d.DAT', datadir, fnums(end)); +else + fn0 = [datadir hIndex.hname{1}(3:end)]; + fn1 = [datadir hIndex.hname{end}(3:end)]; +end +% NOT NEEDED FOR WORKING WITH WAV FILES +% [dtv0,yr,jd,hr,mn,sc,sRate,t,y]=readCF2hdr(fn0,0); +% [dtv1,yr,jd,hr,mn,sc,sRate,t,y]=readCF2hdr(fn1,0); +dtv0 = hIndex.time(1); +dtv1 = hIndex.time(end); +% +if (ischar(baseyear) && strcmp(baseyear, 'auto')) +% baseyear = sub(datevec(dtv0), 1); % new dt encoding from soundIn +baseyear = sub(datevec(hIndex.time(1)), 1); +end +% +% Get start and stop times that include only whole intervals. +dt0 = ceil(dtv0/period) * period; % get start of first whole period +dt1 = floor(dtv1/period) * period; % get end of last whole period +% +% get the system response to subtract out of power vaues later +%[SR] = get_sys_response(PA_rev,gain,filt_sw,sRate,frameSize); +%[SR] = get_sys_response(9,0,4,double(sRate),frameSize); + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% %%%%%%%%%%%%%%%%% Start the loop for calculating LTSpect %%%%%%%%%%%%%% +% +if (1) + gram = zeros((frameSize+zeroPad)/2, round((dt1-dt0)/period)); + fRate = sRate / (frameSize - nOverlap); + printf('%d bins:\n', nCols(gram)) + for i = 1 : nCols(gram) + fprintf(1,'.%s', iff(mod(i,60) == 0, char(10), '')); + dt = dt0 + (i-1) * period; +% Do around 1 million samples at a time. + t0 = dt; + gSum = zeros((frameSize+zeroPad)/2, 1); + ng = 0; +% while (t0 < dt + period) + %t1 = min([dt+period double(t0+1e6/sRate/24/60/60)]); %for big AVG's + t1 = dt + period; + %xt = haruSoundAtTime_CF2_2channel([t0 t1], hIndex, datadir, chan); +% xt = haruSoundAtTime_CF2([t0 t1], hIndex, datadir); %note - calling CF2 specific routine + xt = haruSoundAtTime_sf([t0 t1], hIndex, datadir); + + %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + %xt_2 = xt{1}; %xt_3 = xt_2(:,1); %extra step do to 2 channel data + %x_mn = xt_3 - mean(xt_3); %have to do this in 2 steps now because 2 channel data + x_mn = xt{1} - mean(xt{1}); %remove the sample mean from the chunk of data for 16bit + x = (2.5/(2^16)).*x_mn; % convert to Volts by (Voltrange/bitrange)- for 16bit + %x = (4/(2^12)).*x_mn; %convert to Volts -12bit + g = joespect_LTacoustic(x, frameSize, nOverlap, zeroPad, 'hann', sRate); %5000 for sample rate + +% [g,f,t] = spectrogram(x,hann(frameSize),nOverlap,frameSize,sRate); +% g = abs(g(2:end)); + % last input on line above is sample rate + if (nCols(g) == 0), break; end % skip partial frame at end of day + gSum = gSum + mean(g,2); % sum up the energy - will divide later + ng = ng + nCols(g); + %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + %t0 = t0 + nCols(g)/.9765624 / 24/60/60; %have to change this to make it run m10 hickup - + t0 = t0 + nCols(g) / fRate / 24/60/60; %change back to this + %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% end + grm = 10*log10(gSum); %to convert power to dB + % + %GRAM = grm - SR'; + % + %gram(:,i) = GRAM; %must add this in to get system response removed + gram(:,i) = grm; + % +% + end + fprintf(1, '\n'); +end +% remove the system gain from the spectra +%LTgram = gram - repmat(SR',1,length(gram(1,:))); + +%% +save([datadir 'sg639_GoMex_May18-1kHz-LTSA_1Hz1s.mat'],'gram','dt0','dt1','frameSize','sRate'); + +% \ No newline at end of file diff --git a/dataProcessing/LTSAs/longTermSpectra_NRS.m b/dataProcessing/LTSAs/longTermSpectra_NRS.m new file mode 100644 index 0000000..57d754d --- /dev/null +++ b/dataProcessing/LTSAs/longTermSpectra_NRS.m @@ -0,0 +1,177 @@ +% script modified by J. Haxel May 2015 to calculate the long-term time +% averaged spectra of acoustic signals at NRS +% system response = (preAmp + instResp) is removed and energy values are +% normalized to dB re 1uPa @ 1m +% +% CF2 data - need to use *CF2*.m for all subroutines +% +% data info +%c +if(1) + datadir = '\\ACOUSTIC\data\ONRSN\Samoa\2016-2017\'; + locname = 'NRS10'; + fnums = 1:49999; + chan = 1;a + frameSize = 5000; + fdates = 'C:\Users\haver\Documents\DATAFILE_INDEX\ONRSN\file_dates-NRS10-2016-2017.txt'; +elseif(0) + datadir = '\\ACOUSTIC\data\ONRSN\NRS09\2016-2017\'; + locname = 'NRS09'; + fnums = 1:1750; + chan = 1; + frameSize = 5000; + fdates = 'C:\Users\haver\Documents\DATAFILE_INDEX\ONRSN\file_dates-NRS09-2016-2017.txt'; +elseif(0) + datadir = '\\BIOAC3\datasilo\GLBA\2015\HC-07-Bruise-East\'; + locname = 'HC-07-Bruise-East'; + fnums = 1:3782; % Haru-file numbers + chan = 1; + frameSize = 10000; % for 1Hz resolution @ 5000Hz sample rate + fdates = 'C:\Users\haver\Documents\DATAFILE_INDEX\file_dates-HC-07-Bruise-East-2015.txt'; +elseif(0) + datadir = '\\ACOUSTIC\data\ONRSN\NRS04\2015-2016\'; %'\\ACOUSTIC\data\ONRSN\NRS02\2014-2015'; %'T:\ONRSN\NRS10\2015-2016\' ; + locname = 'NRS04'; + fnums = 1:3136; % Haru-file numbers + chan = 1; + frameSize = 5000; % for 1Hz resolution @ 5000Hz sample rate + fdates = 'C:\Users\haver\Documents\DATAFILE_INDEX\ONRSN\file_dates-NRS04-2015-2016.txt'; +elseif(0) + datadir = '\\ACOUSTIC\data\ONRSN\NRS11\2015-2017\'; %'\\ACOUSTIC\data\ONRSN\NRS02\2014-2015'; %'T:\ONRSN\NRS10\2015-2016\' ; + locname = 'NRS11'; + fnums = 0:4335; % Haru-file numbers + chan = 1; + frameSize = 5000; % for 1Hz resolution @ 5000Hz sample rate + fdates = 'C:\Users\haver\Documents\DATAFILE_INDEX\ONRSN\file_dates-NRS11-2015-2017.txt'; %C:\Users\haver\Documents\DATAFILE_INDEX\ONRSN\file_dates-NRS02-2014-2015.txt'; +elseif(0) + datadir = 'V:\NRS11\2015-2017\'; %V:\NRS11\2015-2017' %'\\ACOUSTIC\data\ONRSN\NRS02\2014-2015'; %'T:\ONRSN\NRS10\2015-2016\' ; + locname = 'NRS11'; + fnums = 0:4335; % Haru-file numbers + chan = 1; + frameSize = 5000; % for 1Hz resolution @ 5000Hz sample rate + fdates = 'C:\Users\haver\Documents\DATAFILE_INDEX\ONRSN\file_dates-NRS11-2015-2017.txt'; %C:\Users\haver\Documents\DATAFILE_INDEX\ONRSN\file_dates-NRS02-2014-2015.txt'; +elseif(0) + datadir = 'T:\ONRSN\NRS06\2014-2015\' %'\\ACOUSTIC\data\ONRSN\NRS02\2014-2015'; %'T:\ONRSN\NRS10\2015-2016\' ; + locname = 'NRS06'; + fnums = 1:3552; % Haru-file numbers + chan = 1; + frameSize = 5000; % for 1Hz resolution @ 5000Hz sample rate + fdates = 'C:\Users\haver\Documents\DATAFILE_INDEX\ONRSN\file_dates-NRS06-2014-2015.txt'; %C:\Users\haver\Documents\DATAFILE_INDEX\ONRSN\file_dates-NRS02-2014-2015.txt'; +elseif(0) + datadir = '\\ACOUSTIC\data\ONRSN\NRS10\2015-2016\'; %'T:\ONRSN\NRS10\2015-2016\' ; + locname = 'NRS10'; + fnums = 0:1775; % Haru-file numbers + chan = 1; + frameSize = 5000; % for 1Hz resolution @ 5000Hz sample rate + fdates = 'C:\Users\haver\Documents\DATAFILE_INDEX\ONRSN\file_dates-NRS10-2015-2016.txt'; +elseif (0) + datadir = '\\ACOUSTIC\data\ANT\RossSea2014\H-11\'; + locname = 'RossSea'; + fnums = 0:739; % Haru-file numbers + frameSize = 1000; % for 1Hz resolution @ 1000Hz sample rate + fdates = 'C:\haxel\DATAFILE_INDEX\ROSS_SEA\file_dates-H-11-RossSea2014.txt'; +elseif (0) + datadir = '\\ACOUSTIC\gloc\Oregon_Coast\SETS_lander\REEF_AprJul2014\'; + locname = 'reef'; + sRate = 32000; % sample rate + fnums = 0:5552; % Haru-file numbers + frameSize = 32000; % for 1Hz resolution @ 32000Hz sample rate + fdates = 'C:\haxel\DATAFILE_INDEX\PMEC_SETS\file_dates-SETSreef-Lander2014.txt'; +elseif (0) + datadir = 'E:\sg639\wav\sg639-1kHz\'; + locname = 'GoMex'; + fnums = 1:49999; + chan = 1; + frameSize = 1000; + fdates = 'C:\Users\haver\Documents\DATAFILE_INDEX\ONRSN\file_dates-NRS10-2016-2017.txt'; +end + % +% Parameters for time averaging and start/stop times +period = 5/1440 %1/24 %5/1440 %hourly =1/24 %600/(60*60*24); % 600 sec in days +baseyear = 'auto'; % auto means get it from first day of data +zeroPad = 0; % zero-pad the signal +nOverlap = 0; % overlap data points +% +%%%%%%%%%%%%%%%%%%%%%%% Configure data Grab %%%%%%%%%%%%%%%%%%%%%%%%%%%% +% +hIndex = readHarufileDateIndex(fdates); +if (0) %change for the old format of datafile + fn0 = sprintf('%sdatafile.%03d', datadir, fnums(1)); + fn1 = sprintf('%sdatafile.%03d', datadir, fnums(end)); +else + fn0 = sprintf('%s%08d.DAT', datadir, fnums(1)); + fn1 = sprintf('%s%08d.DAT', datadir, fnums(end)); +end +[dtv0,yr,jd,hr,mn,sc,sRate,t,y]=readCF2hdr(fn0,0); +[dtv1,yr,jd,hr,mn,sc,sRate,t,y]=readCF2hdr(fn1,0); +% +% +if (ischar(baseyear) && strcmp(baseyear, 'auto')) + baseyear = sub(datevec(dtv0), 1); % new dt encoding from soundIn +end +% +% Get start and stop times that include only whole intervals. +dt0 = ceil(dtv0/period) * period; % get start of first whole period +dt1 = floor(dtv1/period) * period; % get end of last whole period +% +% get the system response to subtract out of power vaues later +%[SR] = get_sys_response(PA_rev,gain,filt_sw,sRate,frameSize); +%[SR] = get_sys_response(9,0,4,double(sRate),frameSize); + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% %%%%%%%%%%%%%%%%% Start the loop for calculating LTSpect %%%%%%%%%%%%%% +% +if (1) + gram = zeros((frameSize+zeroPad)/2, round((dt1-dt0)/period)); + fRate = sRate / (frameSize - nOverlap); + printf('%d dots:', nCols(gram)) + for i = 1 : nCols(gram) + fprintf(1,'.%s', iff(mod(i,60) == 0, char(10), '')); + dt = dt0 + (i-1) * period; +% Do around 1 million samples at a time. + t0 = dt; + gSum = zeros((frameSize+zeroPad)/2, 1); + ng = 0; +% while (t0 < dt + period) + %t1 = min([dt+period double(t0+1e6/sRate/24/60/60)]); %for big AVG's + t1 = dt + period; + %xt = haruSoundAtTime_CF2_2channel([t0 t1], hIndex, datadir, chan); + xt = haruSoundAtTime_CF2([t0 t1], hIndex, datadir); %note - calling CF2 specific routine + %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + %xt_2 = xt{1}; %xt_3 = xt_2(:,1); %extra step do to 2 channel data + %x_mn = xt_3 - mean(xt_3); %have to do this in 2 steps now because 2 channel data + x_mn = xt{1} - mean(xt{1}); %remove the sample mean from the chunk of data for 16bit + x = (2.5/(2^16)).*x_mn; % convert to Volts by (Voltrange/bitrange)- for 16bit + %x = (4/(2^12)).*x_mn; %convert to Volts -12bit + g = joespect_LTacoustic(x, frameSize, nOverlap, zeroPad, 'hann', 5000); %5000 for sample rate + % last input on line above is sample rate + if (nCols(g) == 0), break; end % skip partial frame at end of day + gSum = gSum + mean(g,2); % sum up the energy - will divide later + ng = ng + nCols(g); + %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + %t0 = t0 + nCols(g)/.9765624 / 24/60/60; %have to change this to make it run m10 hickup - + t0 = t0 + nCols(g) / fRate / 24/60/60; %change back to this + %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% end + grm = 10*log10(gSum); %to convert power to dB + % + %GRAM = grm - SR'; + % + %gram(:,i) = GRAM; %must add this in to get system response removed + gram(:,i) = grm; + % +% + end + fprintf(1, '\n'); +end +% remove the system gain from the spectra +%LTgram = gram - repmat(SR',1,length(gram(1,:))); +% +%cd E:\ROSS_SEA-2014 +%cd C:\haxel\OCEAN_NOISE\YAQ_HEAD\ +%cd C:\haxel\LAU_BASIN\ANALYSIS\MATAS +%% +cd C:\Users\haver\Documents\NRS\NRS10\LTSA + +save NRS10_20162017_5minAvgSpect_cal070118 gram dt0 dt1 frameSize sRate + +% \ No newline at end of file diff --git a/dataProcessing/LTSAs/longTermSpectra_NRS_gainAdjust.m b/dataProcessing/LTSAs/longTermSpectra_NRS_gainAdjust.m new file mode 100644 index 0000000..06e32ab --- /dev/null +++ b/dataProcessing/LTSAs/longTermSpectra_NRS_gainAdjust.m @@ -0,0 +1,29 @@ +%% Get NRS system response (SR) from xls pre-amp gain tables +% SR is the system response including hydrophone sensitivity and pre-amp response + + +%read xls preamp info +%[f,g,HS] = importNRScal(workbookFile,sheetName); + +%[f,g,HS] = importNRScal('C:\Users\haver\Documents\NRS\NRS10\NRS10_PreampCal.xlsx','1st Dep'); +%[f,g,HS] = importNRScal_02('C:\Users\haver\Documents\NRS\NRS02\NRS02_PreampCal.xlsx','Dep1'); +[f,g,HS] = importNRScal_2('C:\Users\haver\Documents\NRS\NRS04\NRS04_PreampCal.xlsx','Dep1'); + +%% interpolate values for all Hz +[SR] = calc_sys_response_NRS(f,g,HS); + +%% Apply to gram + +%create matrix of dB adjustment values +dB_gram = repmat(SR, 1, size(gram,2)); + +%% add to gram values and take absolute val +gainAdj_gram = gram - dB_gram; + +%% Check gain! +%gainAdj_gram = gainAdj_gram - 12; + +%% Save LTSA +cd C:\Users\haver\Documents\NRS\NRS04\LTSA + +save NRS04_20152016_1hrAvgSpect_adjusted gainAdj_gram dt0 dt1 frameSize sRate diff --git a/dataProcessing/downsampling/CatBasin2016.m b/dataProcessing/downsampling/CatBasin2016.m new file mode 100644 index 0000000..3cf91e5 --- /dev/null +++ b/dataProcessing/downsampling/CatBasin2016.m @@ -0,0 +1,63 @@ +% convert .flac to .wav and downsample to 10 kHz and 1 kHz + +% need a couple versions of this for different DASBR file formats +% Black DASBRs/SM3M's file format is: +% SM3M-1_0+1_20160717_174310 +% so +% fileOut = [files(f,1).name(14:19) '-' files(f,1).name(21:26) '.wav']; +% + +clear all +clc +warning('off') + +dasbr = true; +glider = false; +tic + +path = 'G:\AFFOGATO-1\DASBR Data - 1\B-4 DASBR\Card-N 512GB\'; +cd(path); +path_out1 = 'E:\BD4-1\BD4-1-256kHz\'; +mkdir(path_out1); +path_out2 = 'E:\BD4-1\BD4-1-10kHz\'; +mkdir(path_out2); +path_out3 = 'E:\BD4-1\BD4-1-1kHz\'; +mkdir(path_out3); + +fs1 = 10000; +fs2 = 1000; + +if (glider) + files = dir('*.flac'); + ch = 1; +else + files = dir('*.wav'); + ch = 2; +end + +for f = 1:length(files) + % fileOut = [files(f,1).name(7:12) '-' files(f,1).name(14:19) '.wav']; % WISPR +% fileOut = [files(f,1).name(14:19) '-' files(f,1).name(21:26) '.wav']; % B1 DASBR + fileOut = [files(f,1).name(13:18) '-' files(f,1).name(20:25) '.wav']; % B2 DASBR + + try + [data,fs] = audioread(files(f,1).name); + + data1 = resample(data,fs1,fs); + data2 = resample(data,fs2,fs); + + audiowrite([path_out1 fileOut],data(:,ch),fs); + audiowrite([path_out2 fileOut],data1(:,ch),fs1); + audiowrite([path_out3 fileOut],data2(:,ch),fs2); + + + % disp([datestr(now) ': ' folder(i,1).name ' file: ' files(j,1).name ' processed']); + catch + disp(['Attention: ' datestr(now) ' file: ' files(f,1).name ' corrupt']); + end + if rem(f,10) == 0; + disp([num2str(f) ' files done']) + end +end + +toc \ No newline at end of file diff --git a/dataProcessing/downsampling/CatBasinHarps.m b/dataProcessing/downsampling/CatBasinHarps.m new file mode 100644 index 0000000..ced74a3 --- /dev/null +++ b/dataProcessing/downsampling/CatBasinHarps.m @@ -0,0 +1,43 @@ +% downsample HARP xwavs from Catalina, but keep as only 30 min files (they +% downsampled but merged into 1 day files and they are too big for ishmael + +clear all +clc +warning('off') + +instr = 'SOCAL_CB_01_02'; +path_in = ['G:\Catalina\data\xwavs\Full bandwidth\' instr '\']; +path_out1 = ['G:\Catalina\data\xwavs\downsampled\' instr '\10kHz\']; +mkdir(path_out1); +path_out2 = ['G:\Catalina\data\xwavs\downsampled\' instr '\1kHz\']; +mkdir(path_out2); + +fs1 = 10000; +fs2 = 1000; + +folderList = dir(path_in); +for fldr = 3:length(folderList) + + fileList = dir([path_in folderList(fldr).name '\*.wav']); + ch = 1; + + for f = 1:length(fileList) + tic + fileOut1 = [fileList(f,1).name(1:end-6) '_10k.wav']; % gets rid of the .x.wav part + fileOut2 = [fileList(f,1).name(1:end-6) '_1k.wav']; % gets rid of the .x.wav part + try + [s,fs] = audioread([path_in folderList(fldr).name '\' fileList(f,1).name]); + + s1 = resample(s,fs1,fs); + s2 = resample(s,fs2,fs); + + audiowrite([path_out1 fileOut1],s1(:,ch),fs1); + audiowrite([path_out2 fileOut2],s2(:,ch),fs2); + catch + disp(['Attention: ' datestr(now) ' file: ' fileList(f,1).name ' corrupt']); + end + clear s s1 s2 + toc + end + fprintf(1,'%s done\n',folderList(fldr).name); +end diff --git a/dataProcessing/downsampling/GoMex2018.m b/dataProcessing/downsampling/GoMex2018.m new file mode 100644 index 0000000..331cf15 --- /dev/null +++ b/dataProcessing/downsampling/GoMex2018.m @@ -0,0 +1,38 @@ +% convert .flac to .wav and downsample to 10 kHz and 1 kHz + +clear all +clc +warning('off') + +path_in ='E:\GoMex2018\flac\'; +% cd(path_in); +path_out1='E:\GoMex2018\wav\sg639-125kHz\'; +mkdir(path_out1); +path_out2='E:\GoMex2018\wav\sg639-10kHz\'; +mkdir(path_out2); +path_out3='E:\GoMex2018\wav\sg639-1kHz\'; +mkdir(path_out3); + +fs1=10000; +fs2=1000; + +files=dir([path_in '*.flac']); +fprintf(1, '%i files:\n',length(files)); +for f=1:length(files) + try + [data,fs] = audioread([path_in files(f,1).name]); + + data1=resample(data,fs1,fs); + data2=resample(data,fs2,fs); + +% audiowrite([path_out1 files(f,1).name(7:end-12) '-' files(f,1).name(14:19) '.wav'],data,fs); + audiowrite([path_out2 files(f,1).name(7:end-12) '-' files(f,1).name(14:19) '.wav'],data1,fs1); + audiowrite([path_out3 files(f,1).name(7:end-12) '-' files(f,1).name(14:19) '.wav'],data2,fs2); + + % disp([datestr(now) ': ' folder(i,1).name ' file: ' files(j,1).name ' processed']); + catch + disp(['Attention: ' datestr(now) ' file: ' files(f,1).name ' corrupt']); + end + fprintf(1,'.') + if rem(f,80) == 0; fprintf(1,'\n%3d ',floor((length(files)-f)/80)); end +end \ No newline at end of file diff --git a/dataProcessing/downsampling/HK462016.m b/dataProcessing/downsampling/HK462016.m new file mode 100644 index 0000000..4759e8f --- /dev/null +++ b/dataProcessing/downsampling/HK462016.m @@ -0,0 +1,9 @@ + + +files=dir('*.flac'); +fs_new=1000; +for i=1:length(files) + [y,fs] = audioread(files(i,1).name); +y_new=resample(y,fs_new,fs); +audiowrite([files(i,1).name(7:12) '-' files(i,1).name(14:19) '.wav'],y_new,fs_new) +end \ No newline at end of file diff --git a/dataProcessing/downsampling/M3R.m b/dataProcessing/downsampling/M3R.m new file mode 100644 index 0000000..d2d86f9 --- /dev/null +++ b/dataProcessing/downsampling/M3R.m @@ -0,0 +1,38 @@ +% convert .flac to .wav and downsample to 10 kHz and 1 kHz +% don't need to convert to wav. going to try to just do it in .flac format + +clear all +clc +warning('off') + +path_in='H:\M3R_randPeriods_downsampled\96kHz_2016\'; +% cd(path); +% path_out1='H:\score\2016\q001-125\'; +% mkdir(path_out1); +% path_out2='H:\score\2016\q001-10\'; +% mkdir(path_out2); +path_out3='H:\M3R_randPeriods_downsampled\1kHz_2016\'; +mkdir(path_out3); + +% fs1=10000; +fs2=1000; + +files=dir([path_in '*.flac']); + +for f=1:length(files) + try + [data,fs] = audioread([path_in files(f,1).name]); + + % data1=resample(data,fs1,fs); + data2=resample(data,fs2,fs); + + % audiowrite([path_out1 files(f,1).name(7:end-12) '-' files(f,1).name(14:19) '.wav'],data,fs); + % audiowrite([path_out2 files(f,1).name(7:end-12) '-' files(f,1).name(14:19) '.wav'],data1,fs1); + audiowrite([path_out3 files(f,1).name(1:7) '01' files(f,1).name(10:end)],data2,fs2); + fprintf(1, '.'); + if (rem(f,80) == 0), fprintf(1, '\n%3d ', floor((length(files)-f)/80)); end + % disp([datestr(now) ': ' folder(i,1).name ' file: ' files(j,1).name ' processed']); + catch + fprintf(1,'\n Attention: %s file: %s corrupt', datestr(now),files(f,1).name); + end +end \ No newline at end of file diff --git a/dataProcessing/downsampling/M3R2015.m b/dataProcessing/downsampling/M3R2015.m new file mode 100644 index 0000000..d6bbad3 --- /dev/null +++ b/dataProcessing/downsampling/M3R2015.m @@ -0,0 +1,38 @@ +% convert .flac to .wav and downsample to 10 kHz and 1 kHz +% don't need to convert to wav. going to try to just do it in .flac format + +clear all +clc +warning('off') + +path_in='D:\score\2015\finTracks\exampleForDave\tk73_M3R_506\fullBandwidth\'; +% cd(path); +% path_out1='H:\score\2016\q001-125\'; +% mkdir(path_out1); +% path_out2='H:\score\2016\q001-10\'; +% mkdir(path_out2); +path_out3='D:\score\2015\finTracks\exampleForDave\tk73_M3R_506\downsampled\'; +mkdir(path_out3); + +% fs1=10000; +fs2=1000; + +files=dir([path_in '*.flac']); + +for f=1:length(files) + try + [data,fs] = audioread([path_in files(f,1).name]); + + % data1=resample(data,fs1,fs); + data2=resample(data,fs2,fs); + + % audiowrite([path_out1 files(f,1).name(7:end-12) '-' files(f,1).name(14:19) '.wav'],data,fs); + % audiowrite([path_out2 files(f,1).name(7:end-12) '-' files(f,1).name(14:19) '.wav'],data1,fs1); + audiowrite([path_out3 files(f,1).name(1:7) '01' files(f,1).name(10:end-5) '.wav'],data2,fs2); + fprintf(1, '.'); + if (rem(f,80) == 0), fprintf(1, '\n%3d ', floor((length(files)-f)/80)); end + % disp([datestr(now) ': ' folder(i,1).name ' file: ' files(j,1).name ' processed']); + catch + fprintf(1,'\n Attention: %s file: %s corrupt', datestr(now),files(f,1).name); + end +end \ No newline at end of file diff --git a/dataProcessing/downsampling/M3R2015_forDE.m b/dataProcessing/downsampling/M3R2015_forDE.m new file mode 100644 index 0000000..a39050e --- /dev/null +++ b/dataProcessing/downsampling/M3R2015_forDE.m @@ -0,0 +1,45 @@ +% convert .flac to .wav and downsample to 10 kHz and 1 kHz +% don't need to convert to wav. going to try to just do it in .flac format + +clear all +clc +warning('off') + +folderList = dir('H:\M3R_flac_2016\'); +for fldr = [4:82] + scorePhone = str2num(folderList(fldr).name); + +path_in=['H:\M3R_flac_2016\' num2str(scorePhone) '\']; +% cd(path); +% path_out1='H:\score\2016\q001-125\'; +% mkdir(path_out1); +% path_out2='H:\score\2016\q001-10\'; +% mkdir(path_out2); +path_out3=['E:\1kHzData\' num2str(scorePhone) '\']; +mkdir(path_out3); +fprintf('Phone %i\n', scorePhone); + +% fs1=10000; +fs2=1000; + +files=dir([path_in '*.flac']); + +for f=1:length(files) + try + [data,fs] = audioread([path_in files(f,1).name]); + + % data1=resample(data,fs1,fs); + data2=resample(data,fs2,fs); + + % audiowrite([path_out1 files(f,1).name(7:end-12) '-' files(f,1).name(14:19) '.wav'],data,fs); + % audiowrite([path_out2 files(f,1).name(7:end-12) '-' files(f,1).name(14:19) '.wav'],data1,fs1); + audiowrite([path_out3 files(f,1).name(1:7) '01' files(f,1).name(10:end-5) '.wav'],data2,fs2); + fprintf(1, '.'); + if (rem(f,80) == 0), fprintf(1, '\n%3d ', floor((length(files)-f)/80)); end + % disp([datestr(now) ': ' folder(i,1).name ' file: ' files(j,1).name ' processed']); + catch + fprintf(1,'\n Attention: %s file: %s corrupt', datestr(now),files(f,1).name); + end +end + +end \ No newline at end of file diff --git a/dataProcessing/downsampling/MIRC.m b/dataProcessing/downsampling/MIRC.m new file mode 100644 index 0000000..5628ce3 --- /dev/null +++ b/dataProcessing/downsampling/MIRC.m @@ -0,0 +1,40 @@ +clear all +clc +warning('off') + +path='H:\GoA_Jul15\sg203\PAM\'; +path_out1='I:\sg203-HF-194kHz\'; +path_out2='I:\sg203-MF-10kHz\'; +path_out3='I:\sg203-LF-1kHz\'; + +fs1=10000; +fs2=1000; + +folder=dir(path); + +for i=3:length(folder) + files=dir([path folder(i,1).name '\*.flac']); + + parfor j=1:length(files) + try + [data,fs,bits] = flacread([path folder(i,1).name '\' files(j,1).name]); + t1=datenum(files(j,1).name(6:end-5),'yymmdd_HHMMSS'); + chk=files(j,1).bytes/length(data); + + if chk>1.5 + disp(['Problem reading file: ' folder(i,1).name '\' files(j,1).name]); + end + + data1=resample(data,fs1,fs); + data2=resample(data,fs2,fs); + + wavwrite(data,fs,bits,[path_out1 files(j,1).name(6:end-12) '-' files(j,1).name(13:end-5) '.wav']); + wavwrite(data1,fs1,bits,[path_out2 files(j,1).name(6:end-12) '-' files(j,1).name(13:end-5) '.wav']); + wavwrite(data2,fs2,bits,[path_out3 files(j,1).name(6:end-12) '-' files(j,1).name(13:end-5) '.wav']); + + disp([datestr(now) ': ' folder(i,1).name ' file: ' files(j,1).name ' processed']); + catch + disp(['Attention: ' datestr(now) ': ' folder(i,1).name ' file: ' files(j,1).name ' corrupt']); + end + end +end \ No newline at end of file diff --git a/dataProcessing/downsampling/SCORE16.m b/dataProcessing/downsampling/SCORE16.m new file mode 100644 index 0000000..df6b85f --- /dev/null +++ b/dataProcessing/downsampling/SCORE16.m @@ -0,0 +1,36 @@ +% convert .flac to .wav and downsample to 10 kHz and 1 kHz + +clear all +clc +warning('off') + +path='E:\sg639\Acoustics\'; +cd(path); +path_out1='HE:\sg639-125\'; +mkdir(path_out1); +path_out2='E:\sg639-10\'; +mkdir(path_out2); +path_out3='E:\sg639-1\'; +mkdir(path_out3); + +fs1=10000; +fs2=1000; + +files=dir('*.flac'); + +for f=1:length(files) + try + [data,fs] = audioread(files(f,1).name); + + data1=resample(data,fs1,fs); + data2=resample(data,fs2,fs); + + audiowrite([path_out1 files(f,1).name(7:end-12) '-' files(f,1).name(14:19) '.wav'],data,fs); + audiowrite([path_out2 files(f,1).name(7:end-12) '-' files(f,1).name(14:19) '.wav'],data1,fs1); + audiowrite([path_out3 files(f,1).name(7:end-12) '-' files(f,1).name(14:19) '.wav'],data2,fs2); + + % disp([datestr(now) ': ' folder(i,1).name ' file: ' files(j,1).name ' processed']); + catch + disp(['Attention: ' datestr(now) ' file: ' files(f,1).name ' corrupt']); + end +end \ No newline at end of file diff --git a/dataProcessing/downsampling/filesizecheck.m b/dataProcessing/downsampling/filesizecheck.m new file mode 100644 index 0000000..3d826c8 --- /dev/null +++ b/dataProcessing/downsampling/filesizecheck.m @@ -0,0 +1,18 @@ +%path_in='H:\score\2016\sg158-LF-1kHz\'; +path_in='H:\score\2016\sg158-MF-10kHz\'; +%path_in='H:\score\2016\sg158-HF-125kHz\'; + +cd(path_in); +files=dir([path_in '\*.wav']); + +data=[]; + +for f=1:650; %length(files); + [y,Fs]=audioread(files(f).name); + data(f,1)=length(y); + data(f,2)=Fs; +end + +for g=2:length(data(:,1)); + data(g-1,3)=data(g,1)-data(g-1,1); +end \ No newline at end of file diff --git a/dataProcessing/downsampling/flac.exe b/dataProcessing/downsampling/flac.exe new file mode 100644 index 0000000..2815fbf Binary files /dev/null and b/dataProcessing/downsampling/flac.exe differ diff --git a/dataProcessing/downsampling/flacread.m b/dataProcessing/downsampling/flacread.m new file mode 100644 index 0000000..a2f45dd --- /dev/null +++ b/dataProcessing/downsampling/flacread.m @@ -0,0 +1,121 @@ +function [Y,FS,NBITS,encoding_info,tag_info] = flacread(FILE,varargin) +%FLACREAD Read FLAC (".flac") sound file. +% Y = OGGREAD(FILE) reads a OGG file specified by the string FILE, +% returning the sampled data in Y. Amplitude values are in the range [-1,+1]. +% +% [Y,FS,NBITS,encoding_info,tag_info] = OGGREAD(FILE) returns the sample rate (FS) in Hertz +% and the number of bits per sample (NBITS) used to encode the +% data in the file. +% +% 'encoding_info' is a string containing information about the mp3 +% encoding used +% +% 'tag_info' is a string containing the tag information of the file +% +% +% Supports two channel or mono encoded data. +% +% See also OGGWRITE, WAVWRITE, AUREAD, AUWRITE. + +persistent cachetest + +a = length(FILE); +if a >= 4 + exten = FILE(a-3:a); + if strcmp(exten,'.flac')==1 + FILE = strcat(FILE,'.flac'); + end +end +if a <= 3 + FILE = strcat(FILE,'.flac'); +end +if exist(FILE,'file') ~= 2 + error(['File not Found: ' FILE]) +end + +%%%%%% Location of the ".exe" Files +s = which('flacread.m'); +ww = findstr('flacread.m',s); +location = s(1:ww-2); + +%%%%Temporary file%%%%%% +tmpdir=''; +tmpdir=tempdir; +[pa,na,ex]=fileparts(FILE); +tmpfile = fullfile(tmpdir,[na ex '.wav']); +cachefile = fullfile(tmpdir,[na ex '.tmp.wav']); + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +%%%%%%%%%%%%%% Info extraction using "ogginfo.exe"%%%%%%%%%%%%%%%% +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +stat_1=0; + +% [stat_1,raw_info] = dos(['"',location,'\ogginfo.exe"',' ', '"',FILE,'"']); +% raw_info_channels_beg = findstr(raw_info,'Channels: '); +% raw_info_channels_end = findstr(raw_info,'Rate: ')-2; +% info_channels = raw_info(raw_info_channels_beg:raw_info_channels_end); +% raw_info_rate_beg = findstr(raw_info,'Rate: '); +% raw_info_rate_end = findstr(raw_info,'Nominal bitrate: ')-1; +% info_rate = strcat(raw_info(raw_info_rate_beg:raw_info_rate_end),' Hz'); +% raw_info_bit_rate = findstr(raw_info,'Nominal bitrate: '); +% raw_info_bit_rate_end = findstr(raw_info,'kb/s'); +% info_bit_rate = raw_info(raw_info_bit_rate+17:raw_info_bit_rate_end-1); +% info_bit_rate = ['Bit Rate: ',num2str(floor(str2num(info_bit_rate))),' Kb/s']; +% encoding_info = {info_channels info_rate info_bit_rate}; +% %%%%% TAG INFO %%%%% +% if isempty(findstr(raw_info,'User comments section follows...')) ~= 1 +% tag_info_beg = findstr(raw_info,'User comments section follows...')+32; +% tag_info_end = findstr(raw_info,'Vorbis stream 1:')-1; +% tag_info = raw_info(tag_info_beg:tag_info_end); +% else +% tag_info = 'No Tag Info'; +% end + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +%%%%%%%%%%%%%% File Decoding using "oggdec.exe" %%%%%%%%%%%%%%%%%% +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +Y=[];FS=NaN;NBITS=NaN; + +if nargin>1 + if isequal(varargin{end},'cache') + %disp(['start "flaccache" /B /min ' location '\flaccache.cmd "' FILE '" "' cachefile '" "' tmpfile '"']); + dos(['start "flaccache" /B /min ' location '\flaccache.cmd "' FILE '" "' cachefile '" "' tmpfile '"']); + cachetest=FILE; + return + end +end + +i=0; +while isequal(cachetest,FILE) && ~exist(tmpfile,'file') + i=i+1; + pause(0.1) + if i==10 + disp(['waiting for cache: ' FILE]) + end + if i==600 + warning(['ignoring cache for:' FILE]) + break + end +end + +if ~exist(tmpfile,'file') + cachetest=[]; + [stat_2,raw_info] = dos(['"' location,'\flac.exe"', ' -d -o "', tmpfile, '" ', '"',FILE,'"']); + if stat_1 || stat_2 + warning(['Problems in flac.exe while decoding file ' FILE]) + end +end + +if exist(tmpfile,'file') + cachetest=[]; + try + [Y,FS,NBITS] = wavread(tmpfile); % Load the data and delete temporary file + catch + try delete(tmpfile); catch warning(['cannot delete ' tmpfile ]); end + error(['Error while reading converted file ' tmpfile ]) + end + try delete(tmpfile); catch warning(['cannot delete ' tmpfile ]); end +else + error(['Error flac.exe could not convert ' FILE ' to .wav']) +end diff --git a/dataProcessing/downsampling/high_pass.m b/dataProcessing/downsampling/high_pass.m new file mode 100644 index 0000000..309f3cc --- /dev/null +++ b/dataProcessing/downsampling/high_pass.m @@ -0,0 +1,26 @@ +function Hd = high_pass +%HIGH_PASS Returns a discrete-time filter object. + +% +% M-File generated by MATLAB(R) 7.8 and the Signal Processing Toolbox 6.11. +% +% Generated on: 07-Jan-2010 10:09:34 +% + +% Butterworth Highpass filter designed using the BUTTER function. + +% All frequency values are in Hz. +Fs = 194162; % Sampling Frequency + +N = 5; % Order +Fc = 5000; % Cutoff Frequency + +% Calculate the zpk values using the BUTTER function. +[z,p,k] = butter(N, Fc/(Fs/2), 'high'); + +% To avoid round-off errors, do not use the transfer function. Instead +% get the zpk representation and convert it to second-order sections. +[sos_var,g] = zp2sos(z, p, k); +Hd = dfilt.df2sos(sos_var, g); + +% [EOF] diff --git a/dataProcessing/downsampling/sg607_SoCal_Feb20.m b/dataProcessing/downsampling/sg607_SoCal_Feb20.m new file mode 100644 index 0000000..208024b --- /dev/null +++ b/dataProcessing/downsampling/sg607_SoCal_Feb20.m @@ -0,0 +1,47 @@ +% convert downsample to 180 kHz 10 kHz and 1 kHz + 5 kHz + +% clear all +% clc +warning('off') + +path='E:\SoCal2020\sg607wav\'; +% cd(path); +% path_out1='E:\SoCal2020\sg607_downsampled\sg607-180kHz\'; +% mkdir(path_out1); +% path_out2='E:\SoCal2020\sg607_downsampled\sg607-10kHz\'; +% mkdir(path_out2); +% path_out3='E:\SoCal2020\sg607_downsampled\sg607-1kHz\'; +% mkdir(path_out3); +path_out4='E:\SoCal2020\sg607_downsampled\sg607-5kHz\'; +mkdir(path_out4); + +% fs1 = 180260; +% fs2 = 10000; +% fs3 = 1000; +fs4 = 5000; + +files = dir([path '*.wav']); +fsTable = []; + +for f = 1:length(files) + try + [data,fs] = audioread([path files(f,1).name]); + fsTable(f,1) = fs; + + % data1 = resample_PMARXL(data, fs1, fs); + % data2 = resample(data, fs2, fs); + % data3 = resample(data, fs3, fs); + data4 = resample(data, fs4, fs); + + % audiowrite([path_out1 files(f,1).name],data,fs1); + % audiowrite([path_out2 files(f,1).name],data2,fs2); + % audiowrite([path_out3 files(f,1).name],data3,fs3); + audiowrite([path_out4 files(f,1).name], data4, fs4); + + fprintf(1, '%s - file #%i: %s processed\n', datestr(now), f, files(f,1).name); + catch + fprintf(1, 'ATTENTION: %s - file #%i: %s corrupt\n', datestr(now), f, files(f,1).name); + end +end + +% save('sampleRatePerFile_sg607.mat','fsTable'); \ No newline at end of file diff --git a/dataProcessing/downsampling/sg639_SoCal_Feb20.m b/dataProcessing/downsampling/sg639_SoCal_Feb20.m new file mode 100644 index 0000000..8e70b8b --- /dev/null +++ b/dataProcessing/downsampling/sg639_SoCal_Feb20.m @@ -0,0 +1,47 @@ +% convert downsample to 180 kHz 10 kHz and 1 kHz + +clear all +% clc +warning('off') + +path='E:\SoCal2020\sg639wav\'; +% cd(path); +% path_out1='F:\SoCal2020\sg639_downsampled\sg639-180kHz\'; +% mkdir(path_out1); +% path_out2='E:\SoCal2020\sg639_downsampled\sg639-10kHz\'; +% mkdir(path_out2); +% path_out3='E:\SoCal2020\sg639_downsampled\sg639-1kHz\'; +% mkdir(path_out3); +path_out4='E:\SoCal2020\sg639_downsampled\sg639-5kHz\'; +mkdir(path_out4); + +% fs1 = 180260; +% fs2 = 10000; +% fs3 = 1000; +fs4 = 5000; + +files=dir([path '*.wav']); +fsTable = []; + +for f=1:length(files) + try + [data,fs] = audioread([path files(f,1).name]); + fsTable(f,1) = fs; + + % data1 = resample_PMARXL(data, fs1, fs); + % data2 = resample(data, fs2, fs); + % data3 = resample(data, fs3, fs); + data4 = resample(data, fs4, fs); + + % audiowrite([path_out1 files(f,1).name],data,fs1); + % audiowrite([path_out2 files(f,1).name],data2,fs2); + % audiowrite([path_out3 files(f,1).name],data3,fs3); + audiowrite([path_out4 files(f,1).name], data4, fs4); + + fprintf(1, '%s - file #%i: %s processed\n', datestr(now), f, files(f,1).name); + catch + fprintf(1, 'ATTENTION: %s - file #%i: %s corrupt\n', datestr(now), f, files(f,1).name); + end +end + +% save('sampleRatePerFile_sg639.mat','fsTable'); \ No newline at end of file diff --git a/dataProcessing/downsampling/work.m b/dataProcessing/downsampling/work.m new file mode 100644 index 0000000..236e372 --- /dev/null +++ b/dataProcessing/downsampling/work.m @@ -0,0 +1,20 @@ +clear all; +clc; + +fin='J:\sg022_HI_194kHz\'; +fout1='J:\sg022_HI_10kHz\'; +fout2='J:\sg022_HI_1kHz\'; + +fsg=dir([fin '*.wav']); + +for i=1:length(fsg); + i + [y,fs,bit]=wavread([fin fsg(i,1).name]); + di=datenum(fsg(i,1).name(7:end-4),'yyyymmdd-HHMMSS'); + d1=resample(y,10000,fs); + d2=resample(y,1000,fs); + wavwrite(d1,10000,16,[fout1 datestr(di,'yymmdd-HHMMSS')]); + wavwrite(d2,1000,16,[fout2 datestr(di,'yymmdd-HHMMSS')]); +end + + \ No newline at end of file diff --git a/dataProcessing/fileConverters/ConvertDAT2WAV/00000020.DAT b/dataProcessing/fileConverters/ConvertDAT2WAV/00000020.DAT new file mode 100644 index 0000000..e69de29 diff --git a/dataProcessing/fileConverters/ConvertDAT2WAV/RWM1-141024-080242.wav b/dataProcessing/fileConverters/ConvertDAT2WAV/RWM1-141024-080242.wav new file mode 100644 index 0000000..18d5776 Binary files /dev/null and b/dataProcessing/fileConverters/ConvertDAT2WAV/RWM1-141024-080242.wav differ diff --git a/dataProcessing/fileConverters/ConvertDAT2WAV/cf2In.m b/dataProcessing/fileConverters/ConvertDAT2WAV/cf2In.m new file mode 100644 index 0000000..5fc82fb --- /dev/null +++ b/dataProcessing/fileConverters/ConvertDAT2WAV/cf2In.m @@ -0,0 +1,107 @@ +function [sams,nChans,sampleSize,sRate,nLeft,dt,endDt] = ... + cf2In(filename, startFrame, nFrames, chans) +%cf2In Read data from Haru's CF2-format file (.dat; from NOAA/PMEL lab). +% +% [sams,nChans,sampleSize,sRate,nLeft,dt] = haruIn(filename, start, n, chans) +% Read Haru Matsumoto's CF2 file format, which is one format of .DAT file. +% Such files typically have names like 00000318.DAT, though there are other +% formats of .DAT files (see haruIn.m). +% +% Input arguments: +% filename name of file +% start starting sample number; default is 0, meaning the +% first sample in the file +% n number of samples to read per channel; default is Inf, +% meaning read the entire file +% chans which channels to read, but since this file format has only +% one channel, this should be 1; channel numbering starts +% at 1; default is 1 +% +% Output arguments: +% sams column vector of samples (or multi-column array, in the +% theoretical case where length(chans) >= 2) +% nChans how many channels are in the file (equals 1 for all +% extant Haru-format files, but who knows about the future) +% sampleSize number of bytes per sample, which is 2 for this format +% sRate sampling rate, Hz +% nLeft number of sample frames remaining in the file after +% reading (use cf2In(filename, 0, 0) to get the total number +% of samples in the file) +% dt date/time (in GMT) that the file starts, encoded as +% in datenum +% endDt date/time (in GMT) of the first sample after the end of the +% file, encoded as in datenum +% +% Note: Reading of .SDAT files is handled by sdatIn.m, which you can call +% directly if you like or do it via this function. +% +% See also haruIn, psmIn, sdatIn. +% +% Dave Mellinger +% Oct. 2015 + +%formerly: [ept,yr,jd,hr,mn,sc,samprate] = readCF2hdr(filename, datayes) + + +%% Documentation from readCF2hdr.m: +% +%[ept,yr,jd,hr,mn,sc,samplerate]=readH3hdr('/dateDir/00000728.DAT',0); %return start time of file. +%[ept,yr,jd,hr,mn,sc,samplerate, t,y]=readH3hdr(''/dateDir/00000728.DAT',1); % return start time and data. +% +% OUTPUT +% file start times from header +% ept is unix time since 1970 +% yr, jd, hr, min, sc of start time +% samprate from header (approx) +% t= dummy time axis for plottting +% y = time series in counts +% +% Note that t and samprate are approximate. +% to get real sample rate difference file start times + +%% Handle the arguments. + +if (nargin < 2), startFrame = 0; end +if (nargin < 3), nFrames = inf; end +if (nargin < 4 || isnan(chans)), chans = 1; end + +if (length(chans) ~= 1 || chans ~= 1) + error('%s: I can handle only one-channel files.', mfilename); +end + +fd = fopen(filename, 'r', 'b'); % open for reading, little-endian fmt +if (fd < 0) + error('Unable to open this sound file for reading:\n%s', filename); +end + +%% Get header stuff. +st1 = fseek(fd, 90, 'bof'); %#ok +timeStr = fread(fd, 64, '*char').'; % GMT as a string +st2 = fseek(fd, 196, 'bof' ); %#ok +sRate = fread(fd, 1, 'long'); % sample rate + +% Decode the GMT string. +yr = str2double(timeStr(1:3)) + 1900; +jd = str2double(timeStr(5:7)); +hr = str2double(timeStr(9:10)); +mn = str2double(timeStr(12:13)); +sc = str2double(timeStr(15:16)) + str2double(timeStr(18:21)) / 1000; +dt = datenum(yr, 1, jd, hr, mn, sc); + +% Other output args. +nChans = 1; % always 1 (so far anyway) for this file format +sampleSize = 2; % always 2 bytes (so far anyway) for this file format +nTotal = (flength(fd) - 256) / sampleSize / nChans;% total sample frames in file +nLeft = nTotal - startFrame; +endDt = dt + nTotal/sRate; + +%% Get the samples. +sams = []; +if (nFrames > 0 && startFrame < nTotal) + st3 = fseek(fd, 256 + startFrame*2, 'bof'); %#ok + sams = fread(fd, nFrames, 'uint16'); +end + + +fclose(fd); + diff --git a/dataProcessing/fileConverters/ConvertDAT2WAV/convertDAT2WAV.m b/dataProcessing/fileConverters/ConvertDAT2WAV/convertDAT2WAV.m new file mode 100644 index 0000000..12409b8 --- /dev/null +++ b/dataProcessing/fileConverters/ConvertDAT2WAV/convertDAT2WAV.m @@ -0,0 +1,36 @@ +clear all; +clc; + +inDir = 'C:\Users\fregosi\Documents\MATLAB\ConvertDAT2WAV\'; +outDir = 'C:\Users\fregosi\Documents\MATLAB\ConvertDAT2WAV\'; + +% Get geo location (East, West, etc.) +geoLoc = 'RWM1'; + +% parpool command is in parallel commputing toolbox. Check if you ahve that +% ver +% look for it, loook at which version. +% parpool started in 2013b, prior to that it was matlabpool +% im running 2013a, so im going to use that. + +% parpool(4) +%matlabpool(4) +%i got admin warnings but I just clicked ok and it went through + +files = dir([inDir '*.DAT']); + +for i = 1 : length(files) +% parfor i = 1 : length(files) + [sams,nChans,sampleSize,sRate,nLeft,dt,endDt] = ... + haruIn([inDir '/' files(i).name]); + % note: endDt seems to be wrong! + sams2=sams-median(sams); + sams3=sams2/(2^16/2); + [yr,mon,day,hr,mn,sec] = datevec(dt); + + outfile = sprintf('%s%s-%02d%02d%02d-%02d%02d%02d.wav', outDir, geoLoc, ... + mod(yr,100), mon, day, hr, mn, round(sec)); + + audiowrite(outfile, sams3, sRate); + disp([datestr(now) ': file ' num2str(i) ' out of ' num2str(length(files)) ' analyzed']); +end diff --git a/dataProcessing/fileConverters/ConvertDAT2WAV/flength.m b/dataProcessing/fileConverters/ConvertDAT2WAV/flength.m new file mode 100644 index 0000000..1af2c13 --- /dev/null +++ b/dataProcessing/fileConverters/ConvertDAT2WAV/flength.m @@ -0,0 +1,23 @@ +function nbytes = flength(filename) +% FLENGTH File length in bytes. +% +% nbytes = flength(filename) +% Return the length in bytes of FILENAME. +% +% nbytes = flength(fid) +% Return the length in bytes of the open file FID. +% +% See also ftell, fseek, fopen, length, size. + +if (isnumeric(filename)) + fid = filename; + k = ftell(fid); + fseek(fid, 0, 'eof'); + nbytes = ftell(fid); + fseek(fid, k, 'bof'); +else + fid = fopen(filename, 'r'); + fseek(fid, 0, 'eof'); + nbytes = ftell(fid); + fclose(fid); +end diff --git a/dataProcessing/fileConverters/ConvertDAT2WAV/haruIn.m b/dataProcessing/fileConverters/ConvertDAT2WAV/haruIn.m new file mode 100644 index 0000000..4b27ded --- /dev/null +++ b/dataProcessing/fileConverters/ConvertDAT2WAV/haruIn.m @@ -0,0 +1,216 @@ +function [sams,nChans,sampleSize,sRate,nLeft,dt,endDt] = ... + haruIn(filename, startFrame, nFrames, chans) +%haruIn Read data from a Haru-format file (.dat, .sdat; from NOAA/PMEL lab). +% +% [sams,nchan,sampleSize,sRate,nLeft,dt] = haruIn(filename, start, n, chans) +% Read Haru Matsumoto's file format (from Chris Fox's lab). Such files +% typically have names like 00000318.DAT or 00000318.SDAT or DATAFILE.318, or +% possibly r0125502.28 or r0125502.28m130. +% +% Input arguments: +% start starting sample number; default is 0, meaning the +% first sample in the file +% n number of samples to read per channel; default is Inf, +% meaning read the entire file +% chans which channels to read; default is all of them; channel +% numbering starts at 1 +% +% Output arguments: +% sams column vector of samples (or multi-column array, if +% chans's length is >= 2) +% nchan how many channels are in the file (equals 1 for all +% extant Haru-format files, but who knows about the future) +% sampleSize number of bytes per sample, 1 or 2 +% sRate sampling rate, Hz +% nLeft number of sample frames remaining in the file after +% reading (use haruIn(filename, 0, 0) to get the total number +% of samples in the file) +% dt date/time (in GMT) that the file starts, encoded as +% in datenum +% endDt date/time (in GMT) of the first sample after the end of the +% file; this output arg is available only for SDAT files +% +% Some Haru-files (ETP before Nov. 99, MAR) do not encode their sampling rate, +% so I have to guess. I'll tell you when I'm guessing. To override this +% guess, do something like this: +% global opHaruSrate; opHaruSrate=112; +% +% Note: Reading of .SDAT files is handled by sdatIn.m, which you can call +% directly if you like or do it via this function. There are also CF2-format +% .DAT files, which are read by cf2In.m. +% +% Dave Mellinger +% 24 Nov 2003 + +% There are two time-date strings in Haru files. The first is from the QTech +% clock, which is generally the more accurate, temperature-corrected one. +% The second is from the Real-Time Clock (RTC), essentially the computer's +% system clock; it's less accurate. +% +% The RTC is set periodically (Matt said maybe once a day) from the QTech +% clock, so the two should never be far apart. However, upon a reboot, the +% QTech clock doesn't know what time it is, so it gets set from the RTC. +% So the QTech clock can have a small quantum jump at that point, equal to +% the amount that the RTC was off. One good way to estimate the jump is +% by looking at the last Haru-file recorded before the reboot, and calculate +% the difference between the two clocks. A better way would probably be to +% estimate the drift rate between the two clocks, then add the appropriate +% amount of drift to that last-file clock difference. I've never done the +% latter. + +global opHaruSrate opHaruSrateWarned + +if (nargin < 2), startFrame = 0; end +if (nargin < 3), nFrames = inf; end +if (nargin < 4), chans = NaN; end + +if (strcmpi(pathExt(filename), 'sdat')) + if (nargin < 4) + [sams,nChans,sampleSize,sRate,nLeft,dt,endDt] = ... + sdatIn(filename, startFrame, nFrames); + else + [sams,nChans,sampleSize,sRate,nLeft,dt,endDt] = ... + sdatIn(filename, startFrame, nFrames, chans); + end + return +end + +fd = fopen(filename, 'r', 'b'); % read-only; big-endian byte order +if (fd < 0) + error(['Can''t open Haru-format file ' filename]); +end + +f12 = char(fread(fd, 12, 'uchar').'); % first 12 bytes +if (strcmpi(f12(1:3), 'BIR')) + shortForm = 0; +else + % Check for CF2 format of Haru file, as indicated by the string 'CF2' at + % byte 152 of the file. + fseek(fd, 152, 'bof'); + cf2str = fread(fd, 2, '*char').'; + if (strcmpi(cf2str, 'cf')) + % CF2-format file. Let cf2In handle it. Don't want to pass fd because + % cf2In has a requirement that the file be opened big-endian. Since that's + % too much to require, cf2In doesn't take fd as an input, just a filename. + fclose(fd); + [sams,nChans,sampleSize,sRate,nLeft,dt,endDt] = ... + cf2In(filename, startFrame, nFrames, chans); + return + end + fseek(fd, 12, 'bof'); + shortForm = 1; + [~, count] = sscanf(f12, 'H%2d%2d%c%3d%c'); + if (count ~= 5) + [~, count] = sscanf(f12, 'h%2d%2d%c%3d%c'); % sometimes it's lowercase + if (count ~= 5) + fclose(fd); + error('Not a Haru-format file.'); + end + end +end + +% Get the date string, parse into dt. +fY = fread(fd, 1, 'uchar'); +% Short-form or regular (non-SDAT) long-form format. +if (shortForm && (fY == 'y' || fY == 'Y')) + % Yet another format for Haru-files, let's call it 'Y format'. + % Example: H0916N043W00Y064:05:26:31:140064:05:26:31:940 + fseek(fd, 13, 'bof'); % where date string starts + dateStr = [f12(11:12) ' ' fread(fd, 16, '*char').']; +elseif (shortForm) + % 'short' format + % example: H1435N043W99 336:17:41:05:530336:17:43:31:920 (at SOF) + % NO, USE FIRST DATE! SECOND ONE IS WRONG. See comment above, or + % mar99\NE\disk3\datafile.249 . + %fscanf(fd, '%d:%d:%d:%d:%3d',5); % skip first date; second is more accurate + %fseek(fd, 30, 'bof'); % where date string starts + k = fread(fd, 16, 'char').'; + dateStr = [f12(11:12) ' ' char(k)]; +else + % 'long' format + fseek(fd, 74, 'bof'); % where date string starts (first one: QTech) + dateStr = char(fread(fd, 42, 'char').'); +end +[dv,count,~,nextIx] = sscanf(dateStr, '%d %d:%d:%d:%d:', 5); +if (count ~= 5) + fclose(fd); + error(['The date string in a Haru-format file is bad: ' dateStr]); +end +% Work around MATLAB bug: Sometimes it doesn't read the last ':'. +if (dateStr(nextIx) == ':'), nextIx = nextIx+1; end +% Only 2 digits of the last 3 are valid, so they're centiseconds. +if (dateStr(nextIx) == '*') + dv(6) = 0; +else + dv(6) = sscanf(dateStr(nextIx : nextIx + 1), '%d') * 10; % *10 to make ms +end +if (dv(1) < 90) % there are no Haru-files from before 1990 + dv(1) = dv(1) + 100; % short-form dates encode year 2000 as '00'. +end +% Make dt from dv, which is [year-1900 yearday hour min sec msec]. +dt = datenum([dv(1)+1900 1 dv(2) dv(3) dv(4) dv(5)+dv(6)/1000]); +%dt = dv; % old encoding: dt was an odd sort of vector + +if (shortForm) + % Short form of header: only date, other parameters are fixed. + nChans = 1; + if (gexist4('opHaruSrate')) + sRate = opHaruSrate; + else + % Atlantic shortForm files were sampled at 110 Hz, ETP ones at 100 Hz. + % I used to think (pre-3/26/2004) Atlantic files were sampled at 112 Hz. + % Examples: + % \\Syspc\hdd3\haruphone\MAR\Mar01_Feb02\16n049w\disk4\datafile.000 + % \\Syspc\hdd1\haruphone\EPR\Nov98\EPR\Nov98-May99\00n095w\disk1\datafile.* + sRate = iff(upper(f12(10)) == 'W' & str2double(f12(7:9)) >= 80, 100, 110); + if (~gexist4('opHaruSrateWarned')), opHaruSrateWarned = 0; end + if (sRate ~= opHaruSrateWarned) + printf(['Warning! Guessing that sample rate of Haru file is %d; ' ... + 'do "help haruIn" to override'], sRate); + opHaruSrateWarned = sRate; + end + end + sampleSize = 1; + precision = 'uchar'; + valueOffset = 128; + headerSize = 46; % bytes +else + % Read rest of long-form (regular or SDAT) header. + %fseek(fd, 106, 'cof'); + fseek(fd, 222, 'bof'); + sams = fread(fd, 3, 'short'); + nChans = sams(1); + sRate = sams(2); + sampleType = sams(3); + if (sampleType ~= 0 && sampleType ~= 2 && sampleType ~= 3) + fclose(fd); + error(['Unknown sample type in Haru-format file ' filename]); + end + sampleSize = iff(sampleType==0, 1, 2); % 0==>1 byte/sam, 2&3==>2 byte/sam + precision = iff(sampleType==0,'uchar',iff(sampleType==2,'short','ushort')); + valueOffset = iff(sampleType==0, 128, iff(sampleType==2, 16384-109, 32768)); + % -109 is right for Albatross bank phone. + headerSize = 256; % bytes +end + +% Determine whether there's a trailing timestamp, and hence length of trailer. +fseek(fd, -17, 'eof'); +[~,count] = fscanf(fd, '%3d:%2d:%2d:%2d:%d'); +trailerLen = iff(count == 5, 46, 0); + +if (nargin < 4), chans = 1:nChans; end + +fseek(fd, 0, 'eof'); +nFramesTotal = (ftell(fd) - headerSize - trailerLen) / nChans / sampleSize; + +% Skip past rest of header. +fseek(fd, headerSize + startFrame * nChans * sampleSize, 'bof'); + +nFrames = min(nFrames, nFramesTotal - startFrame); % prevent reading past EOF +nLeft = nFramesTotal - nFrames - startFrame; + +% Read the data, reshape into sams array. +sams = fread(fd, nFrames * nChans, precision); +fclose(fd); +sams = reshape(sams, nChans, length(sams) / nChans).'; +sams = sams(:,chans) - valueOffset; diff --git a/dataProcessing/fileConverters/ConvertDAT2WAV/pathExt.m b/dataProcessing/fileConverters/ConvertDAT2WAV/pathExt.m new file mode 100644 index 0000000..768a31f --- /dev/null +++ b/dataProcessing/fileConverters/ConvertDAT2WAV/pathExt.m @@ -0,0 +1,15 @@ +function e = pathExt(p) +% ext = pathExt(pathname) +% Given a pathname, return the extension (without its initial '.'). +% If there is no '.' in the filename, return ''. +% +% See also pathDir, pathRoot, pathFile, filesep. + +extchar = '.'; % character that starts an extension + +e = ''; % default result +p = pathFile(p); % get rid of directory part of path +w = find(p == extchar); +if (length(w)) + e = p(w(length(w))+1:length(p)); +end diff --git a/dataProcessing/fileConverters/ConvertDAT2WAV/pathFile.m b/dataProcessing/fileConverters/ConvertDAT2WAV/pathFile.m new file mode 100644 index 0000000..3a03a10 --- /dev/null +++ b/dataProcessing/fileConverters/ConvertDAT2WAV/pathFile.m @@ -0,0 +1,8 @@ +function f = pathFile(p) +% filename = pathFile(pathname) +% Given a pathname, strip off the leading directory names. +% +% See also pathRoot, pathExt, pathDir, pathFileDisk, filesep, fileparts. + +w = [0 find(p == '/' | p == '\')]; +f = p(w(end)+1 : end); diff --git a/dataProcessing/fileConverters/ConvertDAT2WAV/worked.wav b/dataProcessing/fileConverters/ConvertDAT2WAV/worked.wav new file mode 100644 index 0000000..e69de29 diff --git a/dataProcessing/fileConverters/ConvertDAT2WAV/worked2.wav b/dataProcessing/fileConverters/ConvertDAT2WAV/worked2.wav new file mode 100644 index 0000000..e69de29 diff --git a/dataProcessing/fileConverters/QUTR_flacread.m b/dataProcessing/fileConverters/QUTR_flacread.m new file mode 100644 index 0000000..2105ad0 --- /dev/null +++ b/dataProcessing/fileConverters/QUTR_flacread.m @@ -0,0 +1,42 @@ +%converting sg158 flac files + +clear all +clc +glider='SG158'; +lctn='score'; +dplymnt='Dec15 + +path='G:\SeagliderMay2015Newport\'; +path_out1=['G:\SeagliderMay2015Newport\' glider '_wav\']; + +% fs1=10000; +% fs2=1000; + +% folder=dir(path); +folder(1).name=glider; + +for i=1; %3:length(folder) + files=dir([path folder(i,1).name '\*.flac']); + for j=1:length(files) + try + [data,fs,bits] = flacread([path folder(i,1).name '\' files(j,1).name]); + t1=datenum(files(j,1).name(6:end-5),'yymmdd_HHMMSS'); + chk=files(j,1).bytes/length(data); + + if chk>1.5 + disp(['Problem reading file: ' folder(i,1).name '\' files(j,1).name]); + end + +% data1=resample(data,fs1,fs); +% data2=resample(data,fs2,fs); + + wavwrite(data,fs,bits,[path_out1 files(j,1).name(1:end-5) '.wav']); +% wavwrite(data1,fs1,bits,[path_out2 files(j,1).name(1:end-5) '.wav']); +% wavwrite(data2,fs2,bits,[path_out3 files(j,1).name(1:end-5) '.wav']); + + disp([datestr(now) ': ' folder(i,1).name ' file: ' files(j,1).name ' processed']); + catch + disp(['Attention: ' datestr(now) ': ' folder(i,1).name ' file: ' files(j,1).name ' corrupt']); + end + end +end \ No newline at end of file diff --git a/dataProcessing/fileConverters/corruptionCheck.m b/dataProcessing/fileConverters/corruptionCheck.m new file mode 100644 index 0000000..6ca1218 --- /dev/null +++ b/dataProcessing/fileConverters/corruptionCheck.m @@ -0,0 +1,26 @@ +% M3R file corruption check. + +clear all +clc + +% drive 1 - seagate backup 5 TB +% path_in = 'D:\M3R_flac_2015_100sto800s\'; +% path_in = 'E:\M3R_flac_2015_900s\'; +path_in = 'E:\M3R_flac_2016\'; +folderList = dir(path_in); + +for f = 3:length(folderList); + fileList=dir([path_in folderList(f,1).name '\*.flac']); + fprintf(1,'\n%s - %s (%d files):\n', datestr(now), [path_in folderList(f,1).name '\'], length(fileList)); + for g = 1:length(fileList) + try + info = audioinfo([path_in folderList(f,1).name '\' fileList(g,1).name]); +% [s,fs] = audioread([path_in folderList(f,1).name '\' fileList(g,1).name]); + fprintf(1, '.'); + if (rem(g,80) == 0), fprintf(1, '\n%3d ', floor((length(fileList)-g)/80)); end + catch + fprintf(1,'\n Attention: %s file: %s corrupt', datestr(now),fileList(g,1).name); + end + end +% fprintf(1,'%s done',folderList(f,1).name); +end diff --git a/dataProcessing/fileConverters/glider_flacread.m b/dataProcessing/fileConverters/glider_flacread.m new file mode 100644 index 0000000..751b7bd --- /dev/null +++ b/dataProcessing/fileConverters/glider_flacread.m @@ -0,0 +1,49 @@ +%converting glider flac files + +clear all +clc +gldr='sg639'; +lctn='GoMex'; +dplymnt='Jun19'; + +path_flac='E:\GoMex2018\flac\'; +% path_out1=[path gldr '_' lctn '_' dplymnt '_wav\']; +path_out1 = 'E:\GoMex2018\wav\sg639-125kHz\'; +mkdir(path_out1); + +% fs1=10000; +% fs2=1000; + +% folder=dir(path); +folder(1).name=gldr; +n_good=0; +n_corrupt=0; + +for i=1 %3:length(folder) + files=dir([path_flac folder(i,1).name '\*.flac']); +% files = dir([path_flac '\*.flac']); + for j=1:length(files) + try + [data,fs,bits] = audioread([path_flac folder(i,1).name '\' files(j,1).name]); + t1=datenum(files(j,1).name(6:end-5),'_yymmdd_HHMMSS'); + chk=files(j,1).bytes/length(data); + + if chk>1.5 + disp(['Problem reading file: ' folder(i,1).name '\' files(j,1).name]); + end + + % data1=resample(data,fs1,fs); + % data2=resample(data,fs2,fs); + + wavwrite(data,fs,bits,[path_out1 files(j,1).name(1:end-5) '.wav']); + % wavwrite(data1,fs1,bits,[path_out2 files(j,1).name(1:end-5) '.wav']); + % wavwrite(data2,fs2,bits,[path_out3 files(j,1).name(1:end-5) '.wav']); + + disp([datestr(now) ': ' folder(i,1).name ' file: ' files(j,1).name ' processed']); + n_good=n_good+1 + catch + disp(['Attention: ' datestr(now) ': ' folder(i,1).name ' file: ' files(j,1).name ' corrupt']); + n_corrupt=n_corrupt+1 + end + end +end \ No newline at end of file diff --git a/dataProcessing/fileConverters/temp.wav b/dataProcessing/fileConverters/temp.wav new file mode 100644 index 0000000..1be4c56 Binary files /dev/null and b/dataProcessing/fileConverters/temp.wav differ diff --git a/dataProcessing/fileRenamers/copyFiles.m b/dataProcessing/fileRenamers/copyFiles.m new file mode 100644 index 0000000..d6a2e89 --- /dev/null +++ b/dataProcessing/fileRenamers/copyFiles.m @@ -0,0 +1,77 @@ +% copying files for christina to check. based off the csv i made by +% randomizing the 5 min periods + +clear all +clc + +% for 2015 data. +randList = 'C:\Users\selene\OneDrive\projects\AFFOGATO\finWhaleDetectionsComparison\manualMarking\finDetectorPerformanceCheck_2015_20180115.csv'; +copyFileDir = 'E:\M3R_randPeriods_downsampled\'; + +r = readtable(randList); +tic +for f = 1:height(r); + % tic + if ~isnan(r.phone(f)); + randDate = datetime(r.year(f),r.month(f),r.day(f),r.hour(f),r.min(f),0); + randDateEnd = randDate + minutes(5); + fldr = [num2str(r.phone(f)) '\']; + + if r.phone(f) >= 900; + origFileDir = 'E:\M3R_flac_2015_900s\'; + else + origFileDir = 'D:\M3R_flac_2015_100sto800s\'; + end + + fileList = dir([origFileDir fldr '\']); + for g = 3:length(fileList); % start at 3 to skip the two blanks + fileDates(g) = datetime(fileList(g).name(23:37),... + 'InputFormat','yyyyMMdd''T''HHmmss'); + end + idxStart = find(fileDates < randDate,1,'last'); + idxEnd = find(fileDates > randDateEnd,1,'first'); + copyIdxs = idxStart:idxEnd-1; + + for h = copyIdxs + copyfile([origFileDir fldr fileList(h).name],copyFileDir) + end + printf('%i periods done, %i files copied',f,length(copyIdxs)) + else printf('%i is Nan',f) + end + % toc +end +toc + +%% for 2016 data. +randList = 'C:\Users\selene\OneDrive\projects\AFFOGATO\finWhaleDetectionsComparison\manualMarking\finDetectorPerformanceCheck_2016_20180116.csv'; +origFileDir = 'E:\M3R_flac_2016\'; +copyFileDir = 'E:\M3R_randPeriods_downsampled\'; + +r = readtable(randList); +tic +for f = 1:height(r); + % tic + if ~isnan(r.phone(f)); + randDate = datetime(r.year(f),r.month(f),r.day(f),r.hour(f),r.min(f),0); + randDateEnd = randDate + minutes(5); + fldr = [num2str(r.phone(f)) '\']; + + fileList = dir([origFileDir fldr]); + for g = 3:length(fileList); % start at 3 to skip the two blanks + fileDates(g) = datetime(fileList(g).name(23:37),... + 'InputFormat','yyyyMMdd''T''HHmmss'); + end + idxStart = find(fileDates < randDate,1,'last'); + idxEnd = find(fileDates > randDateEnd,1,'first'); + copyIdxs = idxStart:idxEnd-1; + + for h = copyIdxs + copyfile([origFileDir fldr fileList(h).name],copyFileDir) + end + printf('%i periods done, %i files copied',f,length(copyIdxs)) + else printf('%i is Nan',f) + end + % toc +end +toc + diff --git a/dataProcessing/fileRenamers/renameFiles_WISPR.m b/dataProcessing/fileRenamers/renameFiles_WISPR.m new file mode 100644 index 0000000..88eceb0 --- /dev/null +++ b/dataProcessing/fileRenamers/renameFiles_WISPR.m @@ -0,0 +1,27 @@ +%rename WISPR files for 2015/2016 issue +clear all +cd('X:\Gliders\2015_12_21_SCORE_sg158\SG158\') + +%get all the files in the folder +files=dir('*.flac'); + +%loop through each +for f=1:length(files) + % get the file name + oldName=files(f).name; + % example: wispr_150101_000704 + newName=(['wispr_16' oldName(9:end)]); + movefile(oldName,newName); +end +%% Remove WISPR +% I already ran rename_files_dave.m to make dashes instead of _ + +cd('H:\score\2016\sg158_score_Dec15_wav\'); + +files=dir('*.wav'); + +for f=1:length(files) + oldName=files(f).name; + newName=oldName(7:end); + movefile(oldName,newName); +end diff --git a/dataProcessing/fileRenamers/rename_files_dash.m b/dataProcessing/fileRenamers/rename_files_dash.m new file mode 100644 index 0000000..0778d4d --- /dev/null +++ b/dataProcessing/fileRenamers/rename_files_dash.m @@ -0,0 +1,12 @@ +%Convert file names for LTSA in Triton + +files = dir('*.wav'); + +for i=1:length(files) + [pathname,filename,extension] = fileparts(files(i).name); + % the new name + newFilenameA = filename(1:7); + newFilenameB = filename(6:18); + % rename the file + movefile([filename extension], [newFilenameA,'-', newFilenameB,extension]); + end diff --git a/dataProcessing/fileRenamers/rename_files_dash_DASBR.m b/dataProcessing/fileRenamers/rename_files_dash_DASBR.m new file mode 100644 index 0000000..3f72d5b --- /dev/null +++ b/dataProcessing/fileRenamers/rename_files_dash_DASBR.m @@ -0,0 +1,41 @@ +% Convert file names from DASBR recordings for LTSA in Triton +% specific to B Dasbrs - check numbering of old name for W dasbrs) +% path_in='E:\DASBR Data - 1\B-1 DASBR\Card-A 512GB\'; +% path_in='E:\DASBR Data - 1\B-1 DASBR\Card-B 512GB\'; +path_in='J:\B1\'; + +path_out='J:\B1\forLTSA\'; + +files = dir([path_in '*.wav']); +cd(path_in) + +for f=1:length(files) + fprintf(1,'.') + if rem(f,80) == 0; fprintf(1, '\n%3d ',floor((length(files)-f)/80));end + try + oldName = files(f).name; + newName = [path_out oldName(end-18:end-11) '-' oldName(end-9:end)]; + if (~strcmp(oldName,newName)), movefile(oldName, newName); end + + catch + disp(oldName); + end +end + +% check sample rate of all files and move 96kHz files into own folder +cd(path_out) +files = dir([path_out '*.wav']); + +for f = 1:length(files) + try + wavInfo = audioinfo(files(f).name); + if wavInfo.SampleRate == 96000 + movefile(files(f).name,[path_out '96kHzFiles\' files(f).name]); + end + catch + disp(files(f).name) + end +end + + + diff --git a/dataProcessing/fileRenamers/rename_files_dash_SoundTrap.m b/dataProcessing/fileRenamers/rename_files_dash_SoundTrap.m new file mode 100644 index 0000000..0c864b1 --- /dev/null +++ b/dataProcessing/fileRenamers/rename_files_dash_SoundTrap.m @@ -0,0 +1,26 @@ +% Convert file names from DASBR Soundtrap recordings for LTSA in Triton +% file names are 201883689.160720173948.wav + +path_in='D:\AFFOGATO-1\DASBR Data\W-2 DASBR\Soundtrap\'; +path_out=[path_in 'LTSA\']; +[~,~,~] = mkdir(path_out); % the return args prevent an error msg + + +%**sound files are in last folder, name not a date +folder = '201883689\'; + +files = dir([path_in folder '*.wav']); + +for i=1:length(files) + if rem(i,100)==0; + disp(i); + end + try + [y,Fs] = audioread([path_in folder files(i,1).name]); + fname_new=[files(i,1).name(11:16) '-' files(i,1).name(17:end-4) '.wav']; + audiowrite([path_out fname_new],y,Fs); + catch + disp(files(i,1).name); + end +end + diff --git a/dataProcessing/fileRenamers/rename_files_dave.m b/dataProcessing/fileRenamers/rename_files_dave.m new file mode 100644 index 0000000..42747e8 --- /dev/null +++ b/dataProcessing/fileRenamers/rename_files_dave.m @@ -0,0 +1,9 @@ +%RENAMES GLIDER WISPR FILES FOR CREATING LTSA +% LTSA REQUIRES FILE FORMAT: YYMMDD-HHMMSS.wav + +files = dir('*.wav'); +for i = 1:length(files) + nm = files(i).name; + nm(nm == '_') = '-'; + if (~strcmp(files(i).name, nm)), movefile(files(i).name, nm); end +end \ No newline at end of file diff --git a/dataProcessing/fileRenamers/rename_files_gliderDiveProfiles.m b/dataProcessing/fileRenamers/rename_files_gliderDiveProfiles.m new file mode 100644 index 0000000..c17c63f --- /dev/null +++ b/dataProcessing/fileRenamers/rename_files_gliderDiveProfiles.m @@ -0,0 +1,39 @@ +%Rename Glider profile data files for submission to GRIIDC +% 2018 06 25 +% for LADCGEMM project + +%% profile files +path_in = 'T:\GulfOfMexico\2017\forDataSubmission\profiles\'; +cd(path_in) +fileList = dir([path_in '*.nc']); +for i = 1:length(fileList) + orig = fileList(i).name; + new = ['sg607_20170609T182847Z_p' fileList(i).name(5:8) '_delayed.nc']; + if (~strcmp(orig, new)), movefile(orig, new); end +end + +%% acoustic files +% % VERY SLOW +% path_in = 'T:\GulfOfMexico\2017\forDataSubmission\acoustics_netCDF\'; +% cd(path_in) +% fileList = dir([path_in '*.nc']); +% for i = 1:length(fileList) +% orig = fileList(i).name; +% new = [fileList(i).name(1:end-6) '.nc']; +% if length(orig) ~= 25; movefile(orig, new); end +% disp(i) +% end + +% faster +path_in = 'T:\GulfOfMexico\2017\forDataSubmission\acoustics_netCDF\'; +orig_path = 'T:\GulfOfMexico\2017\Acoustics\WAV\'; +origFileList = dir([orig_path '*.wav']); +cd(path_in) +fileList = dir([path_in '*.nc']); +for i = 1:length(fileList) + orig = fileList(i).name; + new = ['sg' fileList(i).name(3:19) origFileList(i).name(12:13) '_delayed.nc']; + [Status, Msg] = FileRename(orig, new); + if Status ~= 1; disp('error'); end +end + diff --git a/fileExchange/contourcs/contourcs.m b/fileExchange/contourcs/contourcs.m new file mode 100644 index 0000000..d8877ea --- /dev/null +++ b/fileExchange/contourcs/contourcs.m @@ -0,0 +1,74 @@ +function Cout = contourcs(varargin) +%CONTOURCS Wrapper to CONTOURS to Obtain Structure Output +% S = CONTOURCS(...) takes the exact same input arguments as the default +% CONTOURC but its output is a struct array with fields: +% +% Level - contour line value +% Length - number of contour line points +% X - X coordinate array of the contour line +% Y - Y coordinate array of the contour line +% +% See also contourc. + +% Version 1.0 (Aug 11, 2010) +% Written by: Takeshi Ikuma +% Revision History: +% - (Aug. 11, 2010) : initial release + +% Run CONTOURC and get output matrix +try + C = contourc(varargin{:}); +catch ME + throwAsCaller(ME); +end + +% Count number of contour segments found (K) +K = 0; +n0 = 1; +while n0<=size(C,2) + K = K + 1; + n0 = n0 + C(2,n0) + 1; +end + +% initialize output struct +el = cell(K,1); +Cout = struct('Level',el,'Length',el,'X',el,'Y',el); + +% fill the output struct +n0 = 1; +for k = 1:K + Cout(k).Level = C(1,n0); + idx = (n0+1):(n0+C(2,n0)); + Cout(k).Length = C(2,n0); + Cout(k).X = C(1,idx); + Cout(k).Y = C(2,idx); + n0 = idx(end) + 1; % next starting index +end + +% Copyright (c)2010, Takeshi Ikuma +% All rights reserved. +% +% Redistribution and use in source and binary forms, with or without +% modification, are permitted provided that the following conditions are +% met: +% +% * Redistributions of source code must retain the above copyright +% notice, this list of conditions and the following disclaimer. +% * Redistributions in binary form must reproduce the above copyright +% notice, this list of conditions and the following disclaimer in the +% documentation and/or other materials provided with the distribution. +% * Neither the names of its contributors may be used to endorse or +% promote products derived from this software without specific prior +% written permission. +% +% THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS +% IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, +% THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR +% PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR +% CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, +% EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, +% PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR +% PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF +% LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING +% NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS +% SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. diff --git a/fileExchange/contourcs/license.txt b/fileExchange/contourcs/license.txt new file mode 100644 index 0000000..427c5ec --- /dev/null +++ b/fileExchange/contourcs/license.txt @@ -0,0 +1,24 @@ +Copyright (c) 2010, Kesh Ikuma +All rights reserved. + +Redistribution and use in source and binary forms, with or without +modification, are permitted provided that the following conditions are +met: + + * Redistributions of source code must retain the above copyright + notice, this list of conditions and the following disclaimer. + * Redistributions in binary form must reproduce the above copyright + notice, this list of conditions and the following disclaimer in + the documentation and/or other materials provided with the distribution + +THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" +AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE +IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE +ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE +LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR +CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF +SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS +INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN +CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) +ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE +POSSIBILITY OF SUCH DAMAGE. diff --git a/fileExchange/inhull/inhull.m b/fileExchange/inhull/inhull.m new file mode 100644 index 0000000..8f077eb --- /dev/null +++ b/fileExchange/inhull/inhull.m @@ -0,0 +1,179 @@ +function in = inhull(testpts,xyz,tess,tol) +% inhull: tests if a set of points are inside a convex hull +% usage: in = inhull(testpts,xyz) +% usage: in = inhull(testpts,xyz,tess) +% usage: in = inhull(testpts,xyz,tess,tol) +% +% arguments: (input) +% testpts - nxp array to test, n data points, in p dimensions +% If you have many points to test, it is most efficient to +% call this function once with the entire set. +% +% xyz - mxp array of vertices of the convex hull, as used by +% convhulln. +% +% tess - tessellation (or triangulation) generated by convhulln +% If tess is left empty or not supplied, then it will be +% generated. +% +% tol - (OPTIONAL) tolerance on the tests for inclusion in the +% convex hull. You can think of tol as the distance a point +% may possibly lie outside the hull, and still be perceived +% as on the surface of the hull. Because of numerical slop +% nothing can ever be done exactly here. I might guess a +% semi-intelligent value of tol to be +% +% tol = 1.e-13*mean(abs(xyz(:))) +% +% In higher dimensions, the numerical issues of floating +% point arithmetic will probably suggest a larger value +% of tol. +% +% DEFAULT: tol = 0 +% +% arguments: (output) +% in - nx1 logical vector +% in(i) == 1 --> the i'th point was inside the convex hull. +% +% Example usage: The first point should be inside, the second out +% +% xy = randn(20,2); +% tess = convhulln(xy); +% testpoints = [ 0 0; 10 10]; +% in = inhull(testpoints,xy,tess) +% +% in = +% 1 +% 0 +% +% A non-zero count of the number of degenerate simplexes in the hull +% will generate a warning (in 4 or more dimensions.) This warning +% may be disabled off with the command: +% +% warning('off','inhull:degeneracy') +% +% See also: convhull, convhulln, delaunay, delaunayn, tsearch, tsearchn +% +% Author: John D'Errico +% e-mail: woodchips@rochester.rr.com +% Release: 3.0 +% Release date: 10/26/06 + +% get array sizes +% m points, p dimensions +p = size(xyz,2); +[n,c] = size(testpts); +if p ~= c + error 'testpts and xyz must have the same number of columns' +end +if p < 2 + error 'Points must lie in at least a 2-d space.' +end + +% was the convex hull supplied? +if (nargin<3) || isempty(tess) + tess = convhulln(xyz); +end +[nt,c] = size(tess); +if c ~= p + error 'tess array is incompatible with a dimension p space' +end + +% was tol supplied? +if (nargin<4) || isempty(tol) + tol = 0; +end + +% build normal vectors +switch p + case 2 + % really simple for 2-d + nrmls = (xyz(tess(:,1),:) - xyz(tess(:,2),:)) * [0 1;-1 0]; + + % Any degenerate edges? + del = sqrt(sum(nrmls.^2,2)); + degenflag = (del<(max(del)*10*eps)); + if sum(degenflag)>0 + warning('inhull:degeneracy',[num2str(sum(degenflag)), ... + ' degenerate edges identified in the convex hull']) + + % we need to delete those degenerate normal vectors + nrmls(degenflag,:) = []; + nt = size(nrmls,1); + end + case 3 + % use vectorized cross product for 3-d + ab = xyz(tess(:,1),:) - xyz(tess(:,2),:); + ac = xyz(tess(:,1),:) - xyz(tess(:,3),:); + nrmls = cross(ab,ac,2); + degenflag = false(nt,1); + otherwise + % slightly more work in higher dimensions, + nrmls = zeros(nt,p); + degenflag = false(nt,1); + for i = 1:nt + % just in case of a degeneracy + % Note that bsxfun COULD be used in this line, but I have chosen to + % not do so to maintain compatibility. This code is still used by + % users of older releases. + % nullsp = null(bsxfun(@minus,xyz(tess(i,2:end),:),xyz(tess(i,1),:)))'; + nullsp = null(xyz(tess(i,2:end),:) - repmat(xyz(tess(i,1),:),p-1,1))'; + if size(nullsp,1)>1 + degenflag(i) = true; + nrmls(i,:) = NaN; + else + nrmls(i,:) = nullsp; + end + end + if sum(degenflag)>0 + warning('inhull:degeneracy',[num2str(sum(degenflag)), ... + ' degenerate simplexes identified in the convex hull']) + + % we need to delete those degenerate normal vectors + nrmls(degenflag,:) = []; + nt = size(nrmls,1); + end +end + +% scale normal vectors to unit length +nrmllen = sqrt(sum(nrmls.^2,2)); +% again, bsxfun COULD be employed here... +% nrmls = bsxfun(@times,nrmls,1./nrmllen); +nrmls = nrmls.*repmat(1./nrmllen,1,p); + +% center point in the hull +center = mean(xyz,1); + +% any point in the plane of each simplex in the convex hull +a = xyz(tess(~degenflag,1),:); + +% ensure the normals are pointing inwards +% this line too could employ bsxfun... +% dp = sum(bsxfun(@minus,center,a).*nrmls,2); +dp = sum((repmat(center,nt,1) - a).*nrmls,2); +k = dp<0; +nrmls(k,:) = -nrmls(k,:); + +% We want to test if: dot((x - a),N) >= 0 +% If so for all faces of the hull, then x is inside +% the hull. Change this to dot(x,N) >= dot(a,N) +aN = sum(nrmls.*a,2); + +% test, be careful in case there are many points +in = false(n,1); + +% if n is too large, we need to worry about the +% dot product grabbing huge chunks of memory. +memblock = 1e6; +blocks = max(1,floor(n/(memblock/nt))); +aNr = repmat(aN,1,length(1:blocks:n)); +for i = 1:blocks + j = i:blocks:n; + if size(aNr,2) ~= length(j), + aNr = repmat(aN,1,length(j)); + end + in(j) = all((nrmls*testpts(j,:)' - aNr) >= -tol,1)'; +end + + + diff --git a/fileExchange/inhull/license.txt b/fileExchange/inhull/license.txt new file mode 100644 index 0000000..80c19d9 --- /dev/null +++ b/fileExchange/inhull/license.txt @@ -0,0 +1,24 @@ +Copyright (c) 2009, John D'Errico +All rights reserved. + +Redistribution and use in source and binary forms, with or without +modification, are permitted provided that the following conditions are +met: + + * Redistributions of source code must retain the above copyright + notice, this list of conditions and the following disclaimer. + * Redistributions in binary form must reproduce the above copyright + notice, this list of conditions and the following disclaimer in + the documentation and/or other materials provided with the distribution + +THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" +AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE +IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE +ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE +LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR +CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF +SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS +INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN +CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) +ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE +POSSIBILITY OF SUCH DAMAGE. diff --git a/fileExchange/utm2deg/license.txt b/fileExchange/utm2deg/license.txt new file mode 100644 index 0000000..91fedee --- /dev/null +++ b/fileExchange/utm2deg/license.txt @@ -0,0 +1,27 @@ +Copyright (c) 2006, Rafael Palacios +All rights reserved. + +Redistribution and use in source and binary forms, with or without +modification, are permitted provided that the following conditions are +met: + + * Redistributions of source code must retain the above copyright + notice, this list of conditions and the following disclaimer. + * Redistributions in binary form must reproduce the above copyright + notice, this list of conditions and the following disclaimer in + the documentation and/or other materials provided with the distribution + * Neither the name of the Univ Pontificia Comillas nor the names + of its contributors may be used to endorse or promote products derived + from this software without specific prior written permission. + +THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" +AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE +IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE +ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE +LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR +CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF +SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS +INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN +CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) +ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE +POSSIBILITY OF SUCH DAMAGE. diff --git a/fileExchange/utm2deg/utm2deg.m b/fileExchange/utm2deg/utm2deg.m new file mode 100644 index 0000000..4e20f01 --- /dev/null +++ b/fileExchange/utm2deg/utm2deg.m @@ -0,0 +1,134 @@ +function [Lat,Lon] = utm2deg(xx,yy,utmzone) +% ------------------------------------------------------------------------- +% [Lat,Lon] = utm2deg(x,y,utmzone) +% +% Description: Function to convert vectors of UTM coordinates into Lat/Lon vectors (WGS84). +% Some code has been extracted from UTMIP.m function by Gabriel Ruiz Martinez. +% +% Inputs: +% x, y , utmzone. +% +% Outputs: +% Lat: Latitude vector. Degrees. +ddd.ddddd WGS84 +% Lon: Longitude vector. Degrees. +ddd.ddddd WGS84 +% +% Example 1: +% x=[ 458731; 407653; 239027; 230253; 343898; 362850]; +% y=[4462881; 5126290; 4163083; 3171843; 4302285; 2772478]; +% utmzone=['30 T'; '32 T'; '11 S'; '28 R'; '15 S'; '51 R']; +% [Lat, Lon]=utm2deg(x,y,utmzone); +% fprintf('%11.6f ',lat) +% 40.315430 46.283902 37.577834 28.645647 38.855552 25.061780 +% fprintf('%11.6f ',lon) +% -3.485713 7.801235 -119.955246 -17.759537 -94.799019 121.640266 +% +% Example 2: If you need Lat/Lon coordinates in Degrees, Minutes and Seconds +% [Lat, Lon]=utm2deg(x,y,utmzone); +% LatDMS=dms2mat(deg2dms(Lat)) +%LatDMS = +% 40.00 18.00 55.55 +% 46.00 17.00 2.01 +% 37.00 34.00 40.17 +% 28.00 38.00 44.33 +% 38.00 51.00 19.96 +% 25.00 3.00 42.41 +% LonDMS=dms2mat(deg2dms(Lon)) +%LonDMS = +% -3.00 29.00 8.61 +% 7.00 48.00 4.40 +% -119.00 57.00 18.93 +% -17.00 45.00 34.33 +% -94.00 47.00 56.47 +% 121.00 38.00 24.96 +% +% Author: +% Rafael Palacios +% Universidad Pontificia Comillas +% Madrid, Spain +% Version: Apr/06, Jun/06, Aug/06 +% Aug/06: corrected m-Lint warnings +%------------------------------------------------------------------------- + +% Argument checking +% +error(nargchk(3, 3, nargin)); %3 arguments required +n1=length(xx); +n2=length(yy); +n3=size(utmzone,1); +if (n1~=n2 || n1~=n3) + error('x,y and utmzone vectors should have the same number or rows'); +end +c=size(utmzone,2); +if (c~=4) + error('utmzone should be a vector of strings like "30 T"'); +end + + + +% Memory pre-allocation +% +Lat=zeros(n1,1); +Lon=zeros(n1,1); + + +% Main Loop +% +for i=1:n1 + if (utmzone(i,4)>'X' || utmzone(i,4)<'C') + fprintf('utm2deg: Warning utmzone should be a vector of strings like "30 T", not "30 t"\n'); + end + if (utmzone(i,4)>'M') + hemis='N'; % Northern hemisphere + else + hemis='S'; + end + + x=xx(i); + y=yy(i); + zone=str2double(utmzone(i,1:2)); + + sa = 6378137.000000 ; sb = 6356752.314245; + +% e = ( ( ( sa ^ 2 ) - ( sb ^ 2 ) ) ^ 0.5 ) / sa; + e2 = ( ( ( sa ^ 2 ) - ( sb ^ 2 ) ) ^ 0.5 ) / sb; + e2cuadrada = e2 ^ 2; + c = ( sa ^ 2 ) / sb; +% alpha = ( sa - sb ) / sa; %f +% ablandamiento = 1 / alpha; % 1/f + + X = x - 500000; + + if hemis == 'S' || hemis == 's' + Y = y - 10000000; + else + Y = y; + end + + S = ( ( zone * 6 ) - 183 ); + lat = Y / ( 6366197.724 * 0.9996 ); + v = ( c / ( ( 1 + ( e2cuadrada * ( cos(lat) ) ^ 2 ) ) ) ^ 0.5 ) * 0.9996; + a = X / v; + a1 = sin( 2 * lat ); + a2 = a1 * ( cos(lat) ) ^ 2; + j2 = lat + ( a1 / 2 ); + j4 = ( ( 3 * j2 ) + a2 ) / 4; + j6 = ( ( 5 * j4 ) + ( a2 * ( cos(lat) ) ^ 2) ) / 3; + alfa = ( 3 / 4 ) * e2cuadrada; + beta = ( 5 / 3 ) * alfa ^ 2; + gama = ( 35 / 27 ) * alfa ^ 3; + Bm = 0.9996 * c * ( lat - alfa * j2 + beta * j4 - gama * j6 ); + b = ( Y - Bm ) / v; + Epsi = ( ( e2cuadrada * a^ 2 ) / 2 ) * ( cos(lat) )^ 2; + Eps = a * ( 1 - ( Epsi / 3 ) ); + nab = ( b * ( 1 - Epsi ) ) + lat; + senoheps = ( exp(Eps) - exp(-Eps) ) / 2; + Delt = atan(senoheps / (cos(nab) ) ); + TaO = atan(cos(Delt) * tan(nab)); + longitude = (Delt *(180 / pi ) ) + S; + latitude = ( lat + ( 1 + e2cuadrada* (cos(lat)^ 2) - ( 3 / 2 ) * e2cuadrada * sin(lat) * cos(lat) * ( TaO - lat ) ) * ( TaO - lat ) ) * ... + (180 / pi); + + Lat(i)=latitude; + Lon(i)=longitude; + +end \ No newline at end of file