smas

ch09_exam_9_11_higher_harmonic_annual_cycle.m

%% Statistical Methods in the Atmospheric Sciences
% Ch09. Time Series
% Example 9.11. A More Complicated Annual Cycle
clc; clear; close all;

%% This section shows how to load data on MATLAB.
% *There are many ways to load data from your raw files, however you may not use traditional codes & loops 
% which can be done by one function. Personally, I recommend you to look at a section of "Data Import and Export"
% in MathWork's manual (https://www.mathworks.com/help/matlab/data-import-and-export.html).
%
% (Please see NOAA Technical Report NWS 41, p.3-4)
% "No measurable precipitation fell during the five-day period centered on the date"
% 1. For each day of the year, for each station, counts were made of the 
%   number of times, out of the 36 years, that 0, 1, centered on the day 
%   in question;
% 2. These overlapping daily totals were smoothed by costructing 9-day
%   running totals;
% 3. The resulting relative frequencies of 0, 1, and of 2  or more days of
%   precipitation in 5, still very noisy, were then subjected to harmonic
%   analysis;
%
% Loading climatological data for the period 1948-1983 at "El Paso"
filename = 'el_paso.csv';                           % GHCND data downloaded from NOAA 
clim = xlsread(filename,'el_paso', 'C2:D13700');    % [date, PRCP]
[yr, mo, day] = datevec(num2str(clim(:,1)), 'yyyymmdd');
delIdx1 = ((mo == 2) & (day == 29)) | (yr == 1947 | yr == 1984);            % Remove leap-days and data in 1948 or 1983
clim(delIdx1,:) = [];                     
tim = clim(:,1);  [yr, mo, day] = datevec(num2str(tim), 'yyyymmdd');
prcp = clim(:,2);                                   % Precipitation
bnPrcp = prcp; bnPrcp(bnPrcp > 0) = 1;              % Binary precipitation (1:rain, 0:no rain)
prcpFreq5 = conv(bnPrcp, ones(1,5), 'same');
runAverDeno = [(ceil(5/2):5), 5*ones(1, length(prcpFreq5) - ceil(5/2)*2), fliplr(ceil(5/2):5)];
prcpFreq5 = prcpFreq5./runAverDeno';
n = length(prcpFreq5);
nyr = n/365;
prcpFreq5Clim = sum(reshape(prcpFreq5, [365, nyr])' == 0, 1)/nyr*100;
prcpFreq5Clim9day = conv(prcpFreq5Clim, ones(1,9), 'same');                 % 9day running mean
runAverDeno = [(ceil(9/2):9), 9*ones(1, length(prcpFreq5Clim) - ceil(9/2)*2), fliplr(ceil(9/2):9)];
prcpFreq5Clim9day = prcpFreq5Clim9day./runAverDeno;


%% Calculating Higher Harmonics for Representing Irregular Time-series
% *I wrote key equations here, but you can see all details in "Chapter 9.4.4. Higher Harmonics" in SMAS.  
% *You can change the number of harmonics (nk) will be used for figure 9.14.
%
% <Harmonic function of periodic time-series>
% yt = y_mean + C1*cos(2pi*t/n-phi_1)
%
% Higher harmonic functions
% yt = y_mean + Sigma_k=1_n/2{C(k)*cos[2pi*kt/n - phi(k)]}
%    = y_mean + Sigma_k=1_n/2{A(k)*cos[2pi*kt/n] + B(k)*sin[2pi*kt/n]}
% 
% When the data values are equally spaced in time with no missing values, 
% the same least-squares parameter values to be obtained more easily and 
% efficiently using
% Ak = 2/n*Sigma_t=1_n{yt*cos(2pi*kt/n)}
% Bk = 2/n*Sigma_t=1_n{yt*sin(2pi*kt/n)}
%
% And,
% Ck = [Ak^2 + Bk^2]^0.5
%        { tan^-1(Bk/Ak),                   Ak > 0
% phik = { tan^-1(Bk/Ak)+-pi , or +- 180,   Ak < 0
%        { pi/2, or 90,                     Ak = 0
%
y = prcpFreq5Clim9day; 
y_mean = mean(y);
t = 1:365;
n = length(t);
nk = 3;                      % Number of harmonics

% Parameter Estimation and Constuction of Harmonic series
A = ones(nk,1); B = A; C = A; phi = A;
har = zeros(nk, length(t));
for k = 1:nk
    % Parameter Estimation
    A(k) = 2/n*sum(y.*cos(2*pi()*k*t./n));
    B(k) = 2/n*sum(y.*sin(2*pi()*k*t./n));
    C(k) = (A(k)^2 + B(k)^2)^0.5;    
    if A(k) > 0
        phi(k) = radtodeg(atan(B(k)/A(k)));
    elseif A(k) < 0
        phi(k) = radtodeg(atan(B(k)/A(k))) + 180;
    else
        phi(k) = 90;
    end
    
    % Construction of Harmonic series (without y_mean)
    har(k,:) = A(k)*cos(2*pi()*k*t./n) + B(k)*sin(2*pi()*k*t./n);
end

% Plotting annual cycle (Figure 9.14)
if 1
    f1 = figure('color', 'w', 'outerposition', [100,100,539,367]);
    p1 = plot(1:365, prcpFreq5Clim9day, '-k'); hold on;
    p2 = plot(1:365, sum(har(:,:),1) + y_mean, '-', 'color', [0 176 80]/255, 'lineWidth', 1);    
    ax1 = gca;
    set(ax1, 'Box', 'on', 'xtick', 0:30:360, 'position', [0.155,0.17,0.79,0.75], 'tickdir', 'in')    
    xlim([0 365])
    ylim([0 100])
    text(230, 93, 'EL PASO, TX', 'fontSize', 10.5)
    text(230, 85, 'ZERO WET DAYS IN 5', 'fontSize', 10.5)
    text(10, 10, sprintf('Sum of the first %d harmonics', nk), 'fontsize', 10,...
        'fontsmoothing', 'on', 'fontangle', 'italic', 'fontweight', 'bold')

    xlabel('DAY OF YEAR', 'fontSize', 10.5)
    ylabel('RELATIVE FREQUENCY [%]', 'fontSize', 10.5)
    grid on
end

% Plotting construction of the smooth curve (Harmonics) (Figure 9.15)
if 1
    f2 = figure('color', 'w', 'outerposition', [100,100,612,485]);
    % (a)
    subplot(2,1,1)
    xAxisMax = 394;
    p1 = plot(1:365, prcpFreq5Clim9day, '-k'); hold on; ax1 = gca;
    plot([0, xAxisMax+1], [0 0], 'k', 'linewidth', 0.15);
    p2 = plot(1:365, har(1,:), '--', 'color', [192 0 0]/255, 'lineWidth', 1);
    p3 = plot(1:365, har(2,:), '-', 'color', [0 112 192]/255, 'lineWidth', 1);
    xTickDay = [75,150,225,300,365];
    xTickDayStr = {'75', '150', '225', '300', '365'};
    xTickRad = [1/2, 1, 3/2, 2]*365/2;
    xTickRadStr = {'\pi/2', '\pi', '3\pi/2', '2\pi'};
    plot(reshape(repmat(xTickDay,[3,1]), length(xTickDay)*3,1),...
        repmat([0,1.2,0]',[length(xTickDay), 1]), '-k', 'linewidth', 0.01)
    plot(reshape(repmat(xTickRad,[3,1]), length(xTickRad)*3,1),...
        repmat([0,-1.2,0]',[length(xTickRad), 1]), '-k', 'linewidth', 0.01)
    t1 = text(xTickDay, ones(1, length(xTickDay))*3.5, xTickDayStr,...
        'FontSize', 11, 'HorizontalAlignment', 'center');    
    t2 = text(xTickRad, ones(1, length(xTickRad))*-3, xTickRadStr,...
        'FontSize', 11, 'HorizontalAlignment', 'center');    
    set(ax1, 'Box', 'off', 'xtick', [], 'XColor',[1 1 1], 'ytick', -15:5:15, 'fontSize', 11,...
        'position', [0.1,0.52,0.86,0.401], 'tickdir', 'out')   
    text(-40, 18, '(a)', 'fontsize', 11)
    annotation('textarrow', [0.47, 0.42], [0.92, 0.87], 'String',...
        sprintf('%.1fcos[2\\pi(2t) / 365 - %.2f]', C(2), degtorad(phi(2))), 'headlength', 4, 'headwidth', 5)
    annotation('textarrow', [0.80, 0.78], [0.54, 0.62], 'String',...
        sprintf('%.1fcos[2\\pit / 365 - %.2f]', C(1), degtorad(phi(1))), 'headlength', 4, 'headwidth', 5)
    xlim([0 xAxisMax])
    ylim([-15.005 15])
    grid on
    
    % (b)
    subplot(2,1,2)
    p4 = plot(1:365, prcpFreq5Clim9day, '-k'); hold on; ax2 = gca;
    p5 = plot(1:365, sum(har(1:2,:),1), '-', 'color', [0 176 80]/255, 'lineWidth', 1);    
    plot([0, xAxisMax+1], [0 0], 'k', 'linewidth', 0.15);plot(reshape(repmat(xTickDay,[3,1]), length(xTickDay)*3,1),...
        repmat([0,1.2,0]',[length(xTickDay), 1]), '-k', 'linewidth', 0.01)
    plot(reshape(repmat(xTickRad,[3,1]), length(xTickRad)*3,1),...
        repmat([0,-1.2,0]',[length(xTickRad), 1]), '-k', 'linewidth', 0.01)
    t3 = text(xTickDay, ones(1, length(xTickDay))*3.5, xTickDayStr,...
        'FontSize', 11, 'HorizontalAlignment', 'center');    
    t4 = text(xTickRad, ones(1, length(xTickRad))*-3, xTickRadStr,...
        'FontSize', 11, 'HorizontalAlignment', 'center');
    set(ax2, 'Box', 'off', 'xtick', [], 'XColor',[1 1 1], 'ytick', -20:10:20, 'fontSize', 11,...
        'position', [0.1,0.047,0.86,0.39], 'tickdir', 'out')   
    text(-40, 23, '(b)', 'fontsize', 11)
    annotation('textarrow', [0.42, 0.48], [0.11, 0.15], 'String',...
        sprintf('%.1fcos[2\\pit / 365 - %.2f] \n +%.1fcos[2\\pi(2t) / 365 - %.2f] ',...
        C(1), degtorad(phi(1)), C(2), degtorad(phi(2))), 'headlength', 4, 'headwidth', 5)
    xlim([0 xAxisMax])
    ylim([-27 23])
    grid on
end

% Printing the Estimated Parameters
fprintf('Estimated Parameters\n')
fprintf('===================================================\n')
fprintf('k\t\t A(k)\t\t B(k)\t\tC(k)\t\tPhi(k)\n')
fprintf('---------------------------------------------------\n')
fprintf('%d\t\t% .3f\t\t% .3f\t\t%.1f%%\t\t%.f\n', [(1:k)', A, B, C, phi]')
fprintf('===================================================\n')