import numpy as np
def deseason_harmonic(dat, K, L):
'''
deseasoned_data, season, beta = deseason_harmonic(dat, K, L)
Subtracts the seasonal cycle (season) from the data (data). Season
is calculated by fitting K harmonics of the annual cycle (as well as the
mean) to the data. Assumes on year has L elements (i.e., 365 for daily data,
73 for pentad data, 52 for weekly data, etc.).
Outputs the deseasonalized data, the season, and the fitting parameters (beta)
Written by Eric Oliver, Dalhousie University, 2007-2011
Adapted from original MATLAB script on 28 November 2012
'''
# ensure that mat is a matrix and a column vector
dat = np.mat(dat)
if dat.shape[1]!=0:
dat = dat.T
# length of time series and generate time vector
n = len(dat)
time = np.mat(np.arange(1,n+1)/(1.*L))
# set up mean and harmonics to fit data
P = np.mat(np.ones((n,1)))
for k in range(1,K+1):
P = np.concatenate((P, np.cos(k*2*np.pi*time.T)), 1)
P = np.concatenate((P, np.sin(k*2*np.pi*time.T)), 1)
# Remove seasonal cycle by harmonic regression
beta = (np.linalg.inv(P.T*P)*P.T)*dat
season = P*beta
dat_ds = dat - season
return dat_ds, season, beta