From b6547c22a07847d8ba42a8da9f6d175c2f418860 Mon Sep 17 00:00:00 2001 From: mkzo <61697292+mkzo@users.noreply.github.com> Date: Fri, 25 Sep 2020 11:53:33 -0400 Subject: [PATCH 1/2] Add files via upload --- LagCorrelation.m | 42 ++++++++++++ groupFourierAnalysis.m | 94 +++++++++++++++++++++++++++ noLagCorrelation.m | 65 +++++++++++++++++++ nonOverlapping.m | 78 ++++++++++++++++++++++ nonOverlappingSummary.m | 86 ++++++++++++++++++++++++ pcaThermistors.m | 140 ++++++++++++++++++++++++++++++++++++++++ scatterBinTrial.m | 18 ++++++ summaryAvgFourier.m | 100 ++++++++++++++++++++++++++++ threeSpikesAnalysis.m | 66 +++++++++++++++++++ tidalAnalysis.m | 35 ++++++++++ 10 files changed, 724 insertions(+) create mode 100644 LagCorrelation.m create mode 100644 groupFourierAnalysis.m create mode 100644 noLagCorrelation.m create mode 100644 nonOverlapping.m create mode 100644 nonOverlappingSummary.m create mode 100644 pcaThermistors.m create mode 100644 scatterBinTrial.m create mode 100644 summaryAvgFourier.m create mode 100644 threeSpikesAnalysis.m create mode 100644 tidalAnalysis.m diff --git a/LagCorrelation.m b/LagCorrelation.m new file mode 100644 index 0000000..0241828 --- /dev/null +++ b/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]); \ No newline at end of file diff --git a/groupFourierAnalysis.m b/groupFourierAnalysis.m new file mode 100644 index 0000000..7ebf4e9 --- /dev/null +++ b/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'); + + + + diff --git a/noLagCorrelation.m b/noLagCorrelation.m new file mode 100644 index 0000000..af85098 --- /dev/null +++ b/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]); diff --git a/nonOverlapping.m b/nonOverlapping.m new file mode 100644 index 0000000..3d8c747 --- /dev/null +++ b/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 + diff --git a/nonOverlappingSummary.m b/nonOverlappingSummary.m new file mode 100644 index 0000000..bba7a48 --- /dev/null +++ b/nonOverlappingSummary.m @@ -0,0 +1,86 @@ +%% 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; + 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; + z(j, :) = 2*abs(X(1:NFFT/2+1)); + %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; + + %z(16:17,2) = 0; %getting rid of outliers to better scales + figure; + heatmap(0:1:24.5,1:19,z,'ColorMap',parula); + caxis([0,1]) + title("FFT Time Series at a 2 week interval, Thermistor 01" + fulltag); + xlabel('Frequency/week'); + ylabel('Power'); + 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 + diff --git a/pcaThermistors.m b/pcaThermistors.m new file mode 100644 index 0000000..f59c812 --- /dev/null +++ b/pcaThermistors.m @@ -0,0 +1,140 @@ +%% Custom Settings +% Range over which to average +avgPeriod = 1946; +columnCluster = false; +% Sampling rate (samples/time period) +Fs = 54486; %freq/week!! +% Range of all thermistors being analyzed +range = [1:24]; +%range(range == 2) = []; +%range(range == 13) = []; +%range(range == 16) = []; +%range(range == 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 +for i = range + tag = 'temperature%02d'; + fulltag = sprintf(tag, i); + %normalizing to get temp. anomaly, reducing 0 frequency bias + data(i,:) = data(i,:) - mean(data(i, :)); + %avg = zeros(1,floor((endSample-startSample)/avgPeriod)); + %for j = 1:floor(((endSample-startSample)/avgPeriod)) + % avg(j) = mean(x(1+(avgPeriod*(j-1)):((avgPeriod*j)))); + %end + + %t = startDate + calendarDuration(0,0,0,0,0,0:11.1:(sampleNumber)*11.1); + %plot(t, x); + %hold on; + + % avg + %avgt = startDate + calendarDuration(0,0,0,0,0,0:avgPeriod*11.1:(sampleNumber)*11.1); + %avgt = avgt(1:length(avg)); + %hold on; + %plot(avgt(1:end),avg(1:end),'-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 + + %xlim([startDate, endDate]); + %xtickformat('dd-MMM-yyyy'); + %title("Sampled Temperature Anomalies of " + fulltag); + %xlabel('Date'); + %ylabel('Temperature Anomalies (C)'); + %legend; + [r, p] = corrcoef(data'); + figure; + heatmap(r); + mycolormap = customcolormap_preset('blue-white-red'); + colormap(mycolormap); + title("Pearson Correlation Coefficients"); + xlabel("Thermistor"); + ylabel("Thermistor"); + caxis([-1 1]); + figure; + heatmap(p); + title("p-values"); + xlabel("Thermistor"); + ylabel("Thermistor"); + caxis([0 1]); + %% PCA + figure; + [coeff, score, latent, ~, explained] = pca(data'); + v = 1:24; + vbls = string(v); + h = biplot(coeff(:,1:3),'VarLabels',vbls); + if (columnCluster) + for k = [6,9,20,22] + h(k).Color = 'y'; % Specify red as the line color + end + for k = [2,13,16,24] + h(k).Color = 'r'; % Specify red as the line color + end + for k = [4,11,18,23] + h(k).Color = 'b'; % Specify red as the line color + end + for k = [1, 14, 15] + h(k).Color = 'k'; % Specify red as the line color + end + for k = [7, 8, 21] + h(k).Color = [0.8500, 0.3250, 0.0980]; % Specify red as the line color + end + for k = [3, 12, 17] + h(k).Color = [0, 0.5, 0]; % Specify red as the line color + end + for k = [5, 10, 19] + h(k).Color = [1 0 1]; % Specify red as the line color + end + else + for k = 1:7 + h(k).Color = 'r'; % Specify red as the line color + end + for k = 8:14 + h(k).Color = 'g'; % Specify red as the line color + end + for k = 15:21 + h(k).Color = 'b'; % Specify red as the line color + end + for k = 22:24 + h(k).Color = 'k'; % Specify red as the line color + end + end + title("Principal Component Analysis of TMPSF Thermistors"); + xlabel("Component 1 ("+explained(1)+"%)"); + ylabel("Component 2 ("+explained(2)+"%)"); + zlabel("Component 3 ("+explained(3)+"%)"); + figure; + hold on; + plot(coeff(:,1),'*-'); + plot(coeff(:,2),'*-'); + plot(coeff(:,3),'*-'); + legend("1","2","3"); + xlabel("Thermistor"); + ylabel("Coefficient"); + title("Principal Component Formulas"); \ No newline at end of file diff --git a/scatterBinTrial.m b/scatterBinTrial.m new file mode 100644 index 0000000..710f568 --- /dev/null +++ b/scatterBinTrial.m @@ -0,0 +1,18 @@ +load("fullData.mat"); +%% REAL CODE +t = t'; +figure; +x = []; +tt = []; +for i = 1:24 + x = [x temp(i,:)]; + tt = [tt t]; +end +h = binscatter(tt, x); +h.NumBins = [200 250]; +ax = gca; +ax.ColorScale = 'log'; +colormap fireice + title("TMPSF Absolute Temperatures 9/29/14-6/26/15"); + xlabel('Date'); + ylabel('Temperatures (C)'); \ No newline at end of file diff --git a/summaryAvgFourier.m b/summaryAvgFourier.m new file mode 100644 index 0000000..3adff6c --- /dev/null +++ b/summaryAvgFourier.m @@ -0,0 +1,100 @@ +%% Custom Settings +% Range over which to average +avgPeriodTable = [1946 3892 7784 15567]; +% Sampling rate (samples/time period) +Fs = 54486; %freq/week!! +% Range of all thermistors being analyzed +range = [1]; +% 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); + dtime = ncread(ncfile,'time'); +end +% fft +i = 1; +dtime = dtime/(60*60*24)+datetime(1900,1,1); %Convert to Matlab time +figure; +for avgPeriod = avgPeriodTable + 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); + avg = zeros(1,floor((endSample-startSample)/avgPeriod)); + for j = 1:floor(((endSample-startSample)/avgPeriod)) + 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); + %figure; + + % actual + plot(dtime, x); + xlim([startDate, endDate]); + xtickformat('dd-MMM-yyyy'); + title("Sampled Temperature Anomalies of " + fulltag); + xlabel('Date'); + ylabel('Temperature Anomalies (C)'); + + % avg + avgt = startDate + calendarDuration(0,0,0,0,0,0:avgPeriod*11.1:(sampleNumber)*11.1); + avgt = avgt(1:length(avg)); + hold on; + %plot(avgt(1:end),avg(1:end),'-s'); + % lines before and after the eruption + %xline(datetime(2014, 9, 29)) %sample 1504547 at Friday, April 10, 2015 + %xline(datetime(2015, 6, 26)); %sample 1964206 at Monday, June 6, 2015 + + % fft graph + %figure; + %plot(f, 2*abs(X(1:NFFT/2+1))); + title("Fourier Frequency Plot of " + fulltag); + xlabel('Period'); + ylabel('Temperature Anomalies (C)'); + ylim([0 0.2]); + + % actual fft + avgFs = Fs/avgPeriod; + avgNFFT=length(avg); %NFFT-point DFT + avgX=(fft(avg,avgNFFT))/avgNFFT; %compute DFT using FFT + avgf = avgFs*(0:(avgNFFT/2))/avgNFFT; + + % avg fft graph + plot(avgf, 2*abs(avgX(1:avgNFFT/2+1))); + hold on; +end + title("Avg Fourier Frequency Plot of " + fulltag); + xlabel('Frequency (1/week)'); + ylabel('Temperature Anomalies (C)'); + axis([0 15 0 0.2]); + legend("6 hours", "12 hours", "24 hours", "48 hours"); + + diff --git a/threeSpikesAnalysis.m b/threeSpikesAnalysis.m new file mode 100644 index 0000000..93bbc4e --- /dev/null +++ b/threeSpikesAnalysis.m @@ -0,0 +1,66 @@ +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) +range = 2; +% preallocation for speed +temp = zeros(24, 2102729); +for i = range + tag = 'temperature%02d'; + fulltag = sprintf(tag, i); + temp(i,:) = ncread(ncfile,fulltag); +end +% fft +for i = range + tag = 'temperature%02d'; + fulltag = sprintf(tag, i); + x = temp(i, 1:752274); + y = temp(i, 752274:1504547); + z = temp(i, 1504547:end); + x = x - mean(x); + y = y - mean(y); + z = z - mean(z); + + NFFT=length(x); %NFFT-point DFT + NFFTy=length(y); + NFFTz=length(z); + X=(fft(x,NFFT))/NFFT; %compute DFT using FFT + Y=(fft(y,NFFTy))/NFFTy; + Z=(fft(z,NFFTy))/NFFTz; + Fs = 54486; %freq/week + fx = Fs*(0:(NFFT/2))/NFFT; + fy = Fs*(0:(NFFTy/2))/NFFTy; + fz = Fs*(0:(NFFTz/2))/NFFTz; + figure; + plot(z); + title("Temperature Anomalies in " + fulltag); + xlabel('Samples'); + ylabel('Temperature (C)'); + hold on; + %xline(1504547); % Friday, April 10, 2015 + %xline(1964206); % Monday, June 6, 2015 + figure; + plot(fx, 2*abs(X(1:NFFT/2+1))); + hold on; + plot(fy, 2*abs(Y(1:NFFTy/2+1))); + plot(fz, 2*abs(Z(1:NFFTz/2+1))); + title("Pre-eruption Power Spectrum, 1st Deployment, Thermistor #"+range); + xlabel('Frequency (1/week)'); + ylabel('Power'); + axis([0 20 0 1]); + legend('09/29/14 - 01/03/15', '01/03/15 - 04/10/15', '04/10/15 - 6/26/15'); + + %X=fftshift((fft(x,NFFT))); + %fVals=(-NFFT/2:NFFT/2-1); %DFT Sample points + + %figure; + %plot(fVals,abs(X)); + %title("Fast Fourier Transform of " + fulltag); %Double Sided FFT - with FFTShift + %xlabel('Wave Number'); + %ylabel('DFT Values'); + %axis([-1 10 0 100]); + +end + + diff --git a/tidalAnalysis.m b/tidalAnalysis.m new file mode 100644 index 0000000..93b4668 --- /dev/null +++ b/tidalAnalysis.m @@ -0,0 +1,35 @@ +%Station 46404 (LLNR 765.4) - WEST ASTORIA - 230 NM West of Astoria, OR +%https://www.ndbc.noaa.gov/station_page.php?station=46404&type=0&startyear=2015&startmonth=06&startday=05&endyear=2015&endmonth=07&endday=18&submit=Submit +file_name = fopen('tidalData.txt', 'r'); +data = fscanf(file_name, '%d %d %d %d %d %d %d %f', inf); +data = reshape(data, 8, 8847)'; +% Data is a 2 dimensional array that stores different tidal values +% in every row, and each row contains the date and water column height +data = data(700:end,:); +water_col = data(:,8); +NFFT=length(water_col); +t = datetime(data(1,1),data(1,2),data(1,3)) + calendarDuration(0,0,0,0,0:15:15*(NFFT-1),0); +avg = mean(water_col); +% water col +%figure; +%plot(t, water_col); +%title("Water Column Heigh of Tidal Buoy 46404"); +%xlabel('Date'); +%ylabel('Water Column Anomalies (m)'); +%axis([0 9000 -30 30]); +water_col = water_col - avg; + +% actual fft +Fs = 672; %number of samples in a week +NFFT=length(water_col); %NFFT-point DFT +X = (fft(water_col))/NFFT; %compute DFT using FFT +f = Fs*(0:(NFFT/2))/NFFT; +% fft graph +plot(f, 2*abs(X(1:NFFT/2+1))); +title("Fourier Frequency Plot of Tidal Buoy 46404"); +xlabel('Frequency (1/week)'); +ylabel('Water Column Anomalies (m)'); +axis([0 30 0 1]); + +fclose(file_name); + From 332c0f14459a451db5bbbe1194431c3bd2a68d1e Mon Sep 17 00:00:00 2001 From: mkzo <61697292+mkzo@users.noreply.github.com> Date: Fri, 2 Oct 2020 20:26:37 -0400 Subject: [PATCH 2/2] Add files via upload --- nonOverlappingSummary.m | 36 ++++++++++++++++-------------------- 1 file changed, 16 insertions(+), 20 deletions(-) diff --git a/nonOverlappingSummary.m b/nonOverlappingSummary.m index bba7a48..b00f8d0 100644 --- a/nonOverlappingSummary.m +++ b/nonOverlappingSummary.m @@ -33,9 +33,9 @@ % fft for i = range clear z; - tag = 'temperature%02d'; + z = zeros(19,54487); + tag = 'Thermistor %02d'; fulltag = sprintf(tag, i); - %period %figure; for j = 1:(length(samples)-1) x = data(i, samples(j):samples(j+1)); @@ -50,7 +50,6 @@ %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}); - hold on; end %title("Fourier Frequency Plot of " + fulltag); %xlabel('Frequency (1/week)'); @@ -59,24 +58,21 @@ %leg = legend(datestr(dates)); %leg.FontSize = 6; - %z(16:17,2) = 0; %getting rid of outliers to better scales + %remove frequencies higher than max_freq) + max_freq = 20; + z = z(:, 1:length(f(f