Add files via upload

LagCorrelation.m

@@ -0,0 +1,42 @@
+load("fullData.mat");
+load("distance.mat");
+startSample = 1;
+[sortedDistance, index] = sort(d(2,:));
+perWeek = 60479/(7*24*60); %samples per period - minute
+groups = floor(30*6/perWeek); %up to sample 1800
+coeff = zeros(4,groups);
+pval = zeros(4, groups);
+j = 1;
+for r = [2 13 16 24] %index
+    for i = 0:groups-1
+        sourceTemp = temp(2,1:perWeek);
+        laggedTemp = temp(r,i*perWeek+1:(i+1)*perWeek);
+        A = [sourceTemp' laggedTemp'];
+        [c, p] = corrcoef(A);
+        coeff(j,i+1) = c(1,2);
+        pval(j,i+1) = p(1,2); 
+    end
+    j = j+1;
+end
+
+figure;
+ylabels = num2cell([2 13 16 24]);
+xlabels = num2cell(startSample:groups);
+h1 = heatmap(xlabels, ylabels, coeff);
+h1.Title = 'TMPSF First Deployment Source #2 - 1 Hour Lag Correlation';
+h1.XLabel = 'Minute';
+h1.YLabel = 'Thermistor';
+mycolormap = customcolormap_preset('blue-white-red');
+colormap(mycolormap);
+caxis([-1 1]);
+
+
+figure;
+ylabels = num2cell([2 13 16 24]);
+xlabels = num2cell(startSample:groups);
+h2 = heatmap(xlabels, ylabels, pval);
+h2.Title = 'TMPSF First Deployment Source #2 - 1 Hour Lag p-values';
+h2.XLabel = 'Minute';
+h2.YLabel = 'Thermistor';
+colormap(parula);
+caxis([0 1]);

groupFourierAnalysis.m

@@ -0,0 +1,94 @@
+%% Custom Settings
+% Range over which to average
+avgPeriod = 1946*4;
+% Sampling rate (samples/time period)
+Fs = 54486; %freq/week!!
+% Range of all thermistors being analyzed
+range = [1 2 13 16 24];
+% Date of first sample being analyzed (entire data set starts at datetime(2014, 9, 29))
+startDate = datetime(2014, 9, 29);
+% Date of last sample being analyzed (entire data set starts at datetime(2015, 7, 1))
+endDate = datetime(2015, 6, 26);
+% Starting and Ending sample (entire data set is 2102729 samples)
+startSample = 1;
+endSample = 2102729;
+% Number of samples being analyzed 
+sampleNumber = endSample-startSample;
+% Useful Sample Numbers and Dates:
+% First = 1 at 9/29/2014
+% Before Eruption = 1504547 at Friday, April 10, 2015
+% After Eruption = 1964206 at Monday, June 6, 2015
+% Last = 2102729 at 6/26/2015
+
+%% Code
+ncfile = 'deployment0001_RS03ASHS-MJ03B-07-TMPSFA301-streamed-tmpsf_sample_20140929T190312-20150626T185957.167762.nc' ; % nc file name
+% To get information about the nc file
+ncinfo(ncfile);
+% to display nc file
+ncdisp(ncfile);
+% preallocation for speed
+data = zeros(24, 2102729);
+for i = range
+    tag = 'temperature%02d';
+    fulltag = sprintf(tag, i);
+    data(i,:) = ncread(ncfile,fulltag);
+end
+% fft
+figure;
+for i = range
+    tag = 'temperature%02d';
+    fulltag = sprintf(tag, i);
+    x = data(i,  startSample:endSample); 
+    %normalizing to get temp. anomaly, reducing 0 frequency bias
+    x = x - mean(x);
+    for j = 1:((endSample-startSample)/avgPeriod)-1
+        avg(j) = mean(x(1+(avgPeriod*(j-1)):((avgPeriod*j))));
+    end
+    % actual fft
+    NFFT=length(x); %NFFT-point DFT      
+    X=(fft(x,NFFT))/NFFT; %compute DFT using FFT 
+    f = Fs*(0:(NFFT/2))/NFFT;
+    t = startDate + calendarDuration(0,0,0,0,0,0:11.1:(sampleNumber)*11.1);
+    subplot(5,1,find(range==i));
+    plot(t, x);
+    xlim([startDate, endDate]);
+    xtickformat('dd-MMM-yyyy');
+    title("Sampled Temperature Anomalies of " + fulltag);
+    xlabel('Date');       
+    ylabel('Anom. (C)');
+    avgt = startDate + calendarDuration(0,0,0,0,0,0:avgPeriod*11.1:(sampleNumber)*11.1);
+    hold on;
+    plot(avgt(1:length(avg)),avg,'-s');
+    % lines before and after the eruption
+    xline(datetime(2015, 4, 10)) %sample 1504547 at Friday, April 10, 2015
+    xline(datetime(2015, 6, 6)); %sample 1964206 at Monday, June 6, 2015
+end
+
+figure;
+for i = range
+    tag = 'temperature%02d';
+    fulltag = sprintf(tag, i);
+    x = data(i,  startSample:endSample); 
+    %normalizing to get temp. anomaly, reducing 0 frequency bias
+    x = x - mean(x);
+    for j = 1:((endSample-startSample)/avgPeriod)-1
+        avg(j) = mean(x(1+(avgPeriod*(j-1)):((avgPeriod*j))));
+    end
+    % actual fft
+    NFFT=length(x); %NFFT-point DFT      
+    X=(fft(x,NFFT))/NFFT; %compute DFT using FFT 
+    f = Fs*(0:(NFFT/2))/NFFT;
+    plot(f, 2*abs(X(1:NFFT/2+1)));
+    hold on;
+end
+    hold on;
+    tidalAnalysis
+    title("Fourier Frequency Plot of 5 Probes");
+    xlabel('Frequency (1/week)');       
+    ylabel('Anom. (C)');
+    axis([0 20 0 0.8]); 
+    legend('1','2','13','16','24','Buoy Data');
+
+
+
+

noLagCorrelation.m

@@ -0,0 +1,65 @@
+% Load data 
+load("fullData.mat");
+load("distance.mat");
+% How to deal with gaps in data? Ignore for larger periods, take into
+% account for shorter periods
+source = 2;
+% Possibly filter out tidal frequencies if they prove too strong
+startTime = t(1);
+endTime = t(end);
+period = datetime(0,0,7); % 1 week
+% Create array of datetimes with period as indicated
+timeLine = [startTime:caldays(7):endTime endTime];
+% Loop through datetime array, calculating corrcoef for each period between
+% each thermistor and the source
+coeff = zeros(24,length(timeLine)-1);
+distance = d(2,:);
+% 8, 9, and 16 have the same distances from 2
+distance(9) = distance(9) + 0.0001;
+distance(16) = distance(16) + 0.0002; %so labeling doesnt return an error
+% 11 and 17 has the same distances from 2
+distance(17) = distance(17) + 0.0001;
+[sortedDistance, index] = sort(distance);
+dates = [];
+%range(source) = [];
+j = 1;
+for i = 1:length(timeLine)-1
+    At = t(logical((t >= timeLine(i)) .* (t < timeLine(i+1))));
+    At.Format = 'dd-MMM-yyyy';
+    dates = [dates (string(At(1)) + "-" + string(At(end)))];
+end
+for r = index
+    for i = 1:length(timeLine)-1
+        At = t(logical((t >= timeLine(i)) .* (t < timeLine(i+1))));
+        startAt = find(t == At(1));
+        endAt = find(t == At(end));
+        laggedTemp = temp(r,startAt:endAt);
+        sourceTemp = temp(source,startAt:endAt);
+        A = [laggedTemp' sourceTemp'];
+        c = corrcoef(A); % columns of A are variables, rows are observations
+        coeff(j,i) = c(1,2);
+    end
+    j = j+1;
+end
+figure;
+ylabels = num2cell(index);
+xlabels = num2cell(1:length(timeLine)-1);
+h1 = heatmap(xlabels, ylabels, coeff);
+h1.Title = 'Correlation Timeline Source #2 - Thermistor and Week';
+h1.XLabel = 'Week';
+h1.YLabel = 'Thermistor';
+mycolormap = customcolormap_preset('blue-white-red');
+colormap(mycolormap);
+caxis([-1 1]);
+
+figure;
+ylabels = num2cell(sortedDistance);
+xlabels = num2cell(dates);
+h2 = heatmap(xlabels, ylabels, coeff);
+h2.Title = 'Correlation Timeline Source #2 - Distances and Date';
+h2.XLabel = 'Date Ranges';
+set(struct(h2).NodeChildren(3), 'XTickLabelRotation', 90);
+h2.YLabel = 'Distances (cm)';
+mycolormap = customcolormap_preset('blue-white-red');
+colormap(mycolormap);
+caxis([-1 1]);

nonOverlapping.m

@@ -0,0 +1,78 @@
+%% Custom Settings
+% Range over which to average (number samples taken in the avg period)
+% Sampling rate (samples/time period)
+Fs = 54486; %freq/week (change requires changing the x axis)
+% Range of all thermistors being analyzed
+range = [1];
+% Periods
+weeks = [2];
+% Number of samples being analyzed 
+sampleNumber = 2102729;
+newcolors = {'#ff0000','#ff4000','#ff8000','#ffbf00','#ffff00','#bfff00','#80ff00', '#40ff00', '#00ff00', '#00ff40', '#00ff80', '#00ffbf', '#00ffff', '#00bfff', '#0080ff', '#0040ff', '#0000ff', '#4000ff', '#8000ff', '#bf00ff', '#ff00ff'};
+
+%% Reading in Data
+    ncfile = 'deployment0001_RS03ASHS-MJ03B-07-TMPSFA301-streamed-tmpsf_sample_20140929T190312-20150626T185957.167762.nc' ; % nc file name
+    % To get information about the nc file
+    ncinfo(ncfile);
+    % to display nc file
+    ncdisp(ncfile);
+    % preallocation for speed
+    data = zeros(24, 2102729);
+    for i = range
+    tag = 'temperature%02d';
+    fulltag = sprintf(tag, i);
+    data(i,:) = ncread(ncfile,fulltag);
+    end
+
+%% Code
+for w = weeks
+    % Sample Array
+    samples = 1:w*54486:2102729; %rounds down number of samples
+    % Dates of Periods Array
+    dates = datetime(2014, 9, 29) + calweeks(1:w:w*length(samples));
+
+    % fft
+    for i = range
+        tag = 'temperature%02d';
+        fulltag = sprintf(tag, i);
+        %period
+        figure; 
+        for j = 1:(length(samples)-1)
+            x = data(i, samples(j):samples(j+1)); 
+            %normalizing to get temp. anomaly, reducing 0 frequency bias
+            x = x - mean(x);
+            % actual fft
+            NFFT=length(x); %NFFT-point DFT      
+            X=(fft(x,NFFT))/NFFT; %compute DFT using FFT 
+            f = Fs*(0:(NFFT/2))/NFFT;
+            % fft graph
+            %str = newcolors(j);
+            %color = sscanf(str(2:end),'%2x%2x%2x',[1 3])/255;
+            plot(f, 2*abs(X(1:NFFT/2+1)), 'Color', newcolors{end-j});
+            hold on;
+        end
+        title("Fourier Frequency Plot of " + fulltag);
+        xlabel('Frequency (1/week)');       
+        ylabel('Power');
+        axis([0 100 0 1]); 
+        leg = legend(datestr(dates));
+        leg.FontSize = 6;
+
+        % actual full plot (all periods)
+        %figure; 
+        %t = dates(1) + calendarDuration(0,0,0,0,0,0:11.1:(2102728)*11.1);
+        %x = data(i, 1:2102729); 
+        %x = x - mean(x);
+        %plot(t, x);
+        %xlim([dates(1), dates(end)]);
+        %xtickformat('dd-MMM-yyyy');
+        %title("Sampled Temperature Anomalies of " + fulltag);
+        %xlabel('Date');       
+        %ylabel('Temperature Anomalies (C)');
+        % lines for each period interval
+        %for d = dates
+        %    xline(datetime(d));
+        %end
+    end
+end
+

nonOverlappingSummary.m

@@ -0,0 +1,82 @@
+%% Custom Settings
+% Range over which to average (number samples taken in the avg period)
+% Sampling rate (samples/time period)
+Fs = 54486; %freq/week (change requires changing the x axis)
+% Range of all thermistors being analyzed
+range = [1,2,13,16,24];
+% Periods
+weeks = [2];
+% Number of samples being analyzed 
+sampleNumber = 2102729;
+newcolors = {'#ff0000','#ff4000','#ff8000','#ffbf00','#ffff00','#bfff00','#80ff00', '#40ff00', '#00ff00', '#00ff40', '#00ff80', '#00ffbf', '#00ffff', '#00bfff', '#0080ff', '#0040ff', '#0000ff', '#4000ff', '#8000ff', '#bf00ff', '#ff00ff'};
+
+%% Reading in Data
+    ncfile = 'deployment0001_RS03ASHS-MJ03B-07-TMPSFA301-streamed-tmpsf_sample_20140929T190312-20150626T185957.167762.nc' ; % nc file name
+    % To get information about the nc file
+    ncinfo(ncfile);
+    % to display nc file
+    ncdisp(ncfile);
+    % preallocation for speed
+    data = zeros(24, 2102729);
+    for i = range
+    tag = 'temperature%02d';
+    fulltag = sprintf(tag, i);
+    data(i,:) = ncread(ncfile,fulltag);
+    end
+
+%% Code
+for w = weeks
+    % Sample Array
+    samples = 1:w*54486:2102729; %rounds down number of samples
+    % Dates of Periods Array
+    dates = datetime(2014, 9, 29) + calweeks(1:w:w*length(samples));
+    % fft
+    for i = range
+        clear z;
+        z = zeros(19,54487);
+        tag = 'Thermistor %02d';
+        fulltag = sprintf(tag, i);
+        %figure; 
+        for j = 1:(length(samples)-1)
+            x = data(i, samples(j):samples(j+1)); 
+            %normalizing to get temp. anomaly, reducing 0 frequency bias
+            x = x - mean(x);
+            % actual fft
+            NFFT=length(x); %NFFT-point DFT      
+            X=(fft(x,NFFT))/NFFT; %compute DFT using FFT 
+            f = Fs*(0:(NFFT/2))/NFFT;
+            % fft graph
+            %str = newcolors(j);
+            %color = sscanf(str(2:end),'%2x%2x%2x',[1 3])/255;
+            z(j, :) = 2*abs(X(1:NFFT/2+1));
+            %plot(f, 2*abs(X(1:NFFT/2+1)), 'Color', newcolors{end-j});
+        end
+        %title("Fourier Frequency Plot of " + fulltag);
+        %xlabel('Frequency (1/week)');       
+        %ylabel('Power');
+        %axis([0 100 0 1]); 
+        %leg = legend(datestr(dates));
+        %leg.FontSize = 6;
+
+        %remove frequencies higher than max_freq)
+        max_freq = 20;
+        z = z(:, 1:length(f(f<max_freq)));
+        figure;
+        h = heatmap(z,'ColorMap', parula);
+        % prctile takes in a vector and calculates the 70th percentile
+        max_color = prctile(reshape(z, [1 numel(z)]),99.5);
+        caxis([0,max_color]);
+        title("Power Spectra at a 2 week interval, " + fulltag);
+        xlabel('Frequency (/week)');  
+        % rounded to nearest tenth
+        x_label = round(f(f<max_freq), 1);
+        h.XData = x_label;
+        ylabel('Time (week of deployment)');
+        %actual full plot (all periods)
+        % lines for each period interval
+        %for d = dates
+        %    xline(datetime(d));
+        %end
+    end
+end
+