WxObserver

Periodograms are computed by zero-extending the IMF modes from 111 to 1024 points, then taking a discrete Fourier transform. The plotted data is the magnitude squared of each DFT value, divided by 1024 and multiplied by 2 (to account for the conjugate-symmetric spectrum). Data is NOT windowed before zero-extending.

Here is the duplication using a log-frequency axis for the periodograms:

duplicated-1

Here is the duplication using a linear period axis (1/frequency):

duplicated-2

Here’s a late addition to the post. I graphed the C3 IMF mode versus sunspot number, just like Willis did in his post. Except I put them on top of each other and scaled the data to be visually similar in amplitude. Also, the IMF is actually the negative of C3 as it seems easier to visually compare that way.

C3-vs-SunspotNumber

I’m not sure if I agree with Willis’ conclusions completely. Yes, the phase between C3 and sunspot number does vary throughout, but there is a 1:1 correlation between extrema in the two curves — every sunspot cycle has a corresponding maxima in the inverted C3 curve.

On the other hand, the amplitudes do NOT track so it would appear that C3 depends on other things besides just sunspots.

I don’t know….it seems like by at least some stretch of the imagination this could be considered a solar signal but there’s clearly a lot more going on there too.

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

WxObserver

California Fires 11/9/2018

Welch’s method applied to Sunspot data

Average Daily Pressure at Hilo (2015)

Tidal Data, Tools and Results Archive

Images