109 lines
5.7 KiB
Python
109 lines
5.7 KiB
Python
def mymatrix_smooth(data_in,smooth_dim=(0,),data_per=(0,),window_len=3,window='step'):
|
|
# from mymatrix_operations import mymatrix_smooth
|
|
import numpy as np
|
|
from mymath import mymath_multiplylist
|
|
# %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
|
|
# general input info %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
|
|
# %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
|
|
shape_data_in = data_in.shape
|
|
n_dim = len(shape_data_in)
|
|
# reshape matrix to get smoothing dimension into 0 and 1 dimension
|
|
rs_beg = np.array(range(0,n_dim))
|
|
rs_end = np.array(range(0,n_dim))
|
|
if n_dim > 1:
|
|
rs_beg[np.array([0,smooth_dim[0]])] = (smooth_dim[0],0)
|
|
if len(smooth_dim) > 1:
|
|
rs_beg[np.array([1,smooth_dim[1]])] = (rs_beg[smooth_dim[1]],rs_beg[1])
|
|
rs_end[np.array([1,smooth_dim[1]])] = (smooth_dim[1],1)
|
|
rs_end[np.array([0,smooth_dim[0]])] = (rs_end[smooth_dim[0]],rs_end[0])
|
|
# transpose data, so that smooth_dim are in the first two dimensions
|
|
data_in = data_in.transpose(rs_beg)
|
|
# %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
|
|
# error checking %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
|
|
# %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
|
|
if len(smooth_dim) > 2:
|
|
exit("smoothing only possible in two dimensions")
|
|
|
|
# %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
|
|
# define window function %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
|
|
# %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
|
|
window_len = int(window_len)
|
|
if window=='step':
|
|
window_num = np.ones(window_len/2+1);
|
|
elif window=='hanning':
|
|
window_num = (np.cos(np.array(range(0,window_len/2+1),dtype=float)/(window_len/2)*np.pi)+1)/2;
|
|
else:
|
|
print 'Defined window function not available for smoothing. No smoothing.'
|
|
window_num = 0;
|
|
|
|
# %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
|
|
# create smoothing matrices for matrix multiplications %%%%%%%%%%%%%%%%%%%%%%
|
|
# %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
|
|
# Allocate smoothing matrix
|
|
if n_dim == 1:
|
|
smooth_ML = 1.
|
|
smooth_MR = np.zeros((shape_data_in[0],shape_data_in[0]),dtype=float)
|
|
np.fill_diagonal(smooth_MR[:], window_num[0])
|
|
else:
|
|
smooth_ML = np.zeros((shape_data_in[smooth_dim[0]],shape_data_in[smooth_dim[0]]),dtype=float)
|
|
if len(smooth_dim) > 1:
|
|
smooth_MR = np.zeros((shape_data_in[smooth_dim[1]],shape_data_in[smooth_dim[1]]),dtype=float)
|
|
np.fill_diagonal(smooth_MR[:], window_num[0])
|
|
else:
|
|
smooth_MR = 1.
|
|
np.fill_diagonal(smooth_ML[:], window_num[0])
|
|
|
|
# Fill smoothing matrix with window_num values
|
|
if n_dim == 1:
|
|
for i_d in range(1,window_len/2+1):
|
|
np.fill_diagonal(smooth_MR[i_d:], window_num[i_d])
|
|
np.fill_diagonal(smooth_MR[:,i_d:], window_num[i_d])
|
|
# modify smooth_matrix in case of periodic data
|
|
if data_per[0] == 1:
|
|
np.fill_diagonal(smooth_MR[(shape_data_in[0]-i_d):], window_num[i_d])
|
|
np.fill_diagonal(smooth_MR[:,(shape_data_in[0]-i_d):], window_num[i_d])
|
|
else:
|
|
for i_d in range(1,window_len/2+1):
|
|
np.fill_diagonal(smooth_ML[i_d:], window_num[i_d])
|
|
np.fill_diagonal(smooth_ML[:,i_d:], window_num[i_d])
|
|
# modify smooth_matrix in case of periodic data
|
|
if data_per[0] == 1:
|
|
np.fill_diagonal(smooth_ML[(shape_data_in[0]-i_d):], window_num[i_d])
|
|
np.fill_diagonal(smooth_ML[:,(shape_data_in[0]-i_d):], window_num[i_d])
|
|
|
|
if len(smooth_dim) > 1:
|
|
for i_d in range(1,window_len/2+1):
|
|
np.fill_diagonal(smooth_MR[i_d:], window_num[i_d])
|
|
np.fill_diagonal(smooth_MR[:,i_d:], window_num[i_d])
|
|
# modify smooth_matrix in case of periodic data
|
|
if data_per[1] == 1:
|
|
np.fill_diagonal(smooth_MR[(shape_data_in[1]-i_d):], window_num[i_d])
|
|
np.fill_diagonal(smooth_MR[:,(shape_data_in[1]-i_d):], window_num[i_d])
|
|
|
|
# %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
|
|
# apply smoothing matrices to input data %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
|
|
# %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
|
|
if n_dim == 1:
|
|
data_in = data_in.dot(smooth_MR)
|
|
data_in = data_in/(np.ones(shape_data_in).dot(smooth_MR))
|
|
elif n_dim == 2:
|
|
data_in = smooth_ML.dot(data_in).dot(smooth_MR)
|
|
data_in = data_in/(smooth_ML.dot(np.ones(shape_data_in)).dot(smooth_MR))
|
|
else:
|
|
shape_data_in_new = (shape_data_in[0],shape_data_in[1],mymath_multiplylist(shape_data_in[2:]))
|
|
data_in = np.reshape(data_in,shape_data_in_new)
|
|
# loop over all other dimensions
|
|
for i_d in range(0,mymath_multiplylist(shape_data_in[2:])):
|
|
# smoothing in first dimension
|
|
if len(smooth_dim):
|
|
data_in[:,:,i_d] = smooth_ML.dot(data_in[:,:,i_d]).dot(smooth_MR)
|
|
data_in[:,:,i_d] = data_in[:,:,i_d]/(smooth_ML.dot(np.ones(shape_data_in[:2])).dot(smooth_MR))
|
|
data_in = np.reshape(data_in,shape_data_in)
|
|
# %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
|
|
# define putput data %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
|
|
# %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
|
|
# transpose data back into original shape
|
|
data_in = data_in.transpose(rs_end)
|
|
|
|
return data_in
|