The tidal data is a CSV file generated by saving Willis’ spreadsheet in CSV format. It is to large to attach to this post but you can create it easily by simply opening Willis’ spreadsheet in Excel and then saving it as a CSV file called “All Tide Data.csv”.
The scripts below were run in Matlab to generate the graphics images, but should also run in freemat. They are attached below in the text, you will need to save each script to a separate .m file for use.
Graphs are shown for 13 different locations which have at least 100 years of uninterrupted tidal data.
Each graph shows two periodograms:
1. The original data is graphed using a solid blue line.
2. The data augmented with fake solar signal is graphed using red circles (no lines between points).
This makes it easier to see where the two sets of data are different.
Notes:
There are graphs for all tidal locations with at least 100 years of uninterrupted data.
I could not find 60 years of uninterrupted data for the location Brest (there was only 56 years there that I could find uninterrupted).
The periodogram used is a modified windowed periodogram using a Hanning window of the same length as the data. Other modifications to the periodogram such as the Welch method were not tried.
fgetl.m
function s=fgetl(fp)
s=”;
gulp=100;
while 1
loc=ftell(fp);
c=fread(fp,gulp,’char’);
cr=find(c == 13);
lf=find(c == 10);
if (numel(cr)+numel(lf)) > 0
last=min([cr;lf])-1;
s=[s,c(1:last)’];
fseek(fp,loc+max([cr(1);lf(1)])+1,’bof’);
break;
end
s=[s,c’];
end
hanning.m
function w=hanning(n)
w=0.5*(1-cos(2*pi*(0:n-1)’/(n-1)));
solarSignal.m
function x=solarSignal(magnitude,periodYears,phaseDegrees,length)
%
% returns a “fake” solar signal with a peak amplitude equal to
% “magnitude”, a period in years of “periodYears” and with an
% initial phase of “phaseDegrees” at the first sample.
% The returned array will contain “length” points with an assumed
% time sample interval of 1/12 year.
%
w=2*pi/(12*periodYears);
phi=phaseDegrees/180*pi;
x=magnitude*cos(w*(0:length-1)’+phi);
x=round(x);
findLongestRuns.m
function [lengths,ranges]=findLongestRuns(tides)
%
% tides is an n-by-m array containing tide data for “n” different
% dates and “m” different locations.
%
% for each location (column) in the tides array, finds the longest
% contiguous string of non-zero data. the length of longest contiguous
%
% data for each location is returned in the m-by-1 lengths array.
%
% the beginning and end indices of the longest data runs are returned
% in the columns of the m-by-2 ranges array.
% ranges(k,1) is the start index of longest run for the k’th location
% ranges(k,2) is the end index for the k’th location
%
% For example, tides(ranges(k,1):ranges(k,2),k) returns the longes
% contiguous run of data for the k’th location.
%
nLocs=size(tides,2);
nDates=size(tides,1);
%
% initialize return value arrays
%
lengths=zeros(nLocs,1);
ranges=zeros(nLocs,2);
for loc=1:nLocs
% find all missing data points (where tides array is zero)
missing=find(tides(:,loc) == 0);
% add points before and after last points to make the computation below
% work.
missing=[0;missing;nDates+1];
% compute the amount of non-zero data between each successive pair of
% missing data points
spans=missing(2:end)-missing(1:end-1);
% from this, calculate the maximum un-interrupted span of data and get
% its beginning and ending indices
[lengths(loc),idx]=max(spans);
ranges(loc,:)=[missing(idx)+1,missing(idx+1)-1];
end
loadAllTideData.m
function [year,locations,tides]=loadAllTideData(csvFileName)
%
% loads tidal data from CSV file
%
% year is n-by-1 array containing the year and fractional year of each entry
% locations is a m-by-1 cell array containing strings from the column headers
% tides is an n-by-m array containing tide data for each location columns.
% missing tide data is represented by zeros
a=csvread(csvFileName,1,0);
year=a(:,1);
tides=a(:,2:end);
%
% read header line as a string
fp=fopen(csvFileName,’r’);
s=fgetl(fp);
fclose(fp);
%
% extract station names by breaking at comma characters
%
idx=find(s==’,’)’;
idx=[idx;numel(s)+1];
locations=cell(numel(idx)-1,0);
for k=2:numel(idx)
locations{k-1,1}=s(idx(k-1)+1:idx(k)-1);
end
fakeSignalDemo.m
% this script loads tidal data if not already loaded, then
% generates two FFT’s of the data for one of the locations
% in the tidal data array (determined by variable “loc”).
%
% A “fake” solar signal with period of “solarPeriod” years is added
% to the tidal data before generating the second FFT.
freemat=exist(‘periodogram’,’file’) == 0;
if ~exist(‘year’,’var’)
[year,locations,tides]=loadAllTideData(‘All Tide Data.csv’);
end
if ~exist(‘contiguousLengths’,’var’)
[contiguousLengths,contiguousRanges]=findLongestRuns(tides);
end
requiredLengthYears=100;
goodLocs=find(contiguousLengths >= (requiredLengthYears*12));
%
% loop through each tidal location which has an adequately long
% string of uninterrupted data, and plot two periodograms for that
% location:
% 1) Periodogram of the raw tide data
% 2) Periodogram of raw data plus the fake solar signal
%
% The raw data is plotted with a blue solid line and the augmented
% periodogram is plotted with red circles on top of the blue line
% this makes it easy to see places where the two sets of data are
% different.
%
for locIdx=1:numel(goodLocs)
loc=goodLocs(locIdx); % set this to the index of a tide gauge location to be analyzed
% extract the data subset for this location
x=tides(contiguousRanges(loc,1):contiguousRanges(loc,2),loc);
n=numel(x);
% create time axis in years starting at zero; there are 12 points per
% year.
t=(0:n-1)’/12;
win=hanning(n);
%
% before computing the periodogram the following steps are taken:
% 1) Remove any linear trend by performing a straight line fit and
% subtracting the resulting fit from the data.
% 2) Offset the data so that it has zero mean.
px=polyfit(t,x,1);
x1=x-polyval(px,t);
x2=x1-mean(x1);
%
% create the fake solar signal with amplitude of 2.5mm as specified
% in Dr. Shaviv’s published article. Period in years and phase offset
% may be adjusted to examine effects.
%
solarMag=2.5; % mm
solarPeriod=11.5; % years
solarPhase=0.0; % degrees
solar=solarSignal(solarMag,solarPeriod,solarPhase,n);
y=x+solar;
py=polyfit(t,y,1);
y1=y-polyval(py,t);
y2=y1-mean(y1);
%
% compute the periodogram
%
if freemat
% freemat does not have the periodogram function
% so we synthesize if from an FFT instead
U=win’*win;
Xmag=abs(fft(x2)).^2/(n*U);
Ymag=abs(fft(x2)).^2/(n*U);
omega=2*pi*(0:n-1)’/n;
m=floor(n/2)+1;
Xmag=Xmag(1:m);
Ymag=Ymag(1:m);
omega=omega(1:m);
else
[Xmag,omega]=periodogram(x2,win,n);
[Ymag,omega]=periodogram(y2,win,n);
end
%
% remove zero frequency data point so when we compute period,
% there is no divide by zero error
%
Xmag=Xmag(2:end);
Ymag=Ymag(2:end);
omega=omega(2:end);
%
% compute frequency in years from radian frequency based on
% a sample rate of 12 times per year.
%
freq=omega/2/pi*12;
% scale periodograms to be relative to peak value in percent
Xmag=Xmag/max(Xmag)*100;
Ymag=Ymag/max(Ymag)*100;
% compute period in years
period=1./freq;
%
% plot the data
%
figure(1);
semilogx(period,Xmag,’b’, period,Ymag,’ro’);
axis([0.4,30,0,104]);
grid on;
xlabel(‘Period, years’);
ylabel(‘Amplitude, percent of peak’);
legend(‘Original Data’,’Data plus 11-year Signal’);
title(sprintf(‘Periodogram for %1.0f years of Tide Data at %s’,…
n/12,locations{loc}));
drawnow;
pngfn=sprintf(‘%s.png’,locations{loc});
if freemat
% can print to file in freemat too but not sure what cmd format
% is…
% this pause feature may not work so well in freemat.
% consider removing it and just plotting one location at a time
% by putting a break here or removing the containing loop
pause(3);
else
% uncomment next line to save graph images
print(‘-dpng’,’-r300′,pngfn);
pause;
end
end
