diff --git a/pywafo/src/wafo/covariance/core.py b/pywafo/src/wafo/covariance/core.py index 5a917f7..cc69d9c 100644 --- a/pywafo/src/wafo/covariance/core.py +++ b/pywafo/src/wafo/covariance/core.py @@ -130,8 +130,8 @@ class CovData1D(WafoData): CovData """ - def __init__(self,*args,**kwds): - super(CovData1D, self).__init__(*args,**kwds) + def __init__(self, *args, **kwds): + super(CovData1D, self).__init__(*args, **kwds) self.name = 'WAFO Covariance Object' self.type = 'time' @@ -143,7 +143,7 @@ class CovData1D(WafoData): self.norm = 0 somekeys = ['phi', 'name', 'h', 'tr', 'lagtype', 'v', 'type', 'norm'] - self.__dict__.update(sub_dict_select(kwds,somekeys)) + self.__dict__.update(sub_dict_select(kwds, somekeys)) self.setlabels() def setlabels(self): @@ -153,10 +153,10 @@ class CovData1D(WafoData): ''' N = len(self.type) - if N==0: + if N == 0: raise ValueError('Object does not appear to be initialized, it is empty!') - labels = ['','ACF',''] + labels = ['', 'ACF', ''] if self.lagtype.startswith('t'): labels[0] = 'Lag [s]' @@ -252,36 +252,36 @@ class CovData1D(WafoData): if rate is None: rate = 1 ##interpolation rate else: - rate = 2**nextpow2(rate) ##make sure rate is a power of 2 + rate = 2 ** nextpow2(rate) ##make sure rate is a power of 2 ## add a nugget effect to ensure that round off errors ## do not result in negative spectral estimates - ACF[0] = ACF[0] +nugget + ACF[0] = ACF[0] + nugget n = ACF.size # embedding a circulant vector and Fourier transform if fast: - nfft = 2**nextpow2(2*n-2) + nfft = 2 ** nextpow2(2 * n - 2) else: - nfft = 2*n-2 + nfft = 2 * n - 2 - nf = nfft/2 ## number of frequencies - ACF = r_[ACF,zeros(nfft-2*n+2),ACF[n-1:0:-1]] + nf = nfft / 2 ## number of frequencies + ACF = r_[ACF, zeros(nfft - 2 * n + 2), ACF[n - 1:0:-1]] - Rper = (fft(ACF,nfft).real).clip(0) ## periodogram + Rper = (fft(ACF, nfft).real).clip(0) ## periodogram RperMax = Rper.max() - Rper = where(Rper 1: - So.args = linspace(0, pi/dT, nf*rate) - if method=='stineman': + So.args = linspace(0, pi / dT, nf * rate) + if method == 'stineman': So.data = stineman_interp(So.args, w, S) else: intfun = interpolate.interp1d(w, S, kind=method) @@ -300,11 +300,11 @@ class CovData1D(WafoData): [s] if lagtype=='t' [m] otherwise ''' - dt1 = self.args[1]-self.args[0] - n = size(self.args)-1 - t = self.args[-1]-self.args[0] - dt = t/n - if abs(dt-dt1) > 1e-10: + dt1 = self.args[1] - self.args[0] + n = size(self.args) - 1 + t = self.args[-1] - self.args[0] + dt = t / n + if abs(dt - dt1) > 1e-10: warnings.warn('Data is not uniformly sampled!') return dt @@ -385,7 +385,7 @@ class CovData1D(WafoData): dT = self.sampling_period() - x = zeros((ns, cases+1)) + x = zeros((ns, cases + 1)) if derivative: xder = x.copy() @@ -399,22 +399,22 @@ class CovData1D(WafoData): ## Covariance matrix floatinfo = finfo(float) if (abs(ACF[-1]) > floatinfo.eps): ## assuming ACF(n+1)==0 - m2 = 2*n-1 - nfft = 2**nextpow2(max(m2, 2*ns)) - ACF = r_[ACF, zeros((nfft-m2,1)), ACF[-1:0:-1,:]] + m2 = 2 * n - 1 + nfft = 2 ** nextpow2(max(m2, 2 * ns)) + ACF = r_[ACF, zeros((nfft - m2, 1)), ACF[-1:0:-1, :]] #disp('Warning: I am now assuming that ACF(k)=0 ') #disp('for k>MAXLAG.') else: # # ACF(n)==0 - m2 = 2*n-2 - nfft = 2**nextpow2(max(m2, 2*ns)) - ACF = r_[ACF, zeros((nfft-m2, 1)), ACF[n-1:1:-1, :]] + m2 = 2 * n - 2 + nfft = 2 ** nextpow2(max(m2, 2 * ns)) + ACF = r_[ACF, zeros((nfft - m2, 1)), ACF[n - 1:1:-1, :]] ##m2=2*n-2 - S = fft(ACF,nfft,axis=0).real ## periodogram + S = fft(ACF, nfft, axis=0).real ## periodogram I = S.argmax() - k = flatnonzero(S<0) - if k.size>0: + k = flatnonzero(S < 0) + if k.size > 0: #disp('Warning: Not able to construct a nonnegative circulant ') #disp('vector from the ACF. Apply the parzen windowfunction ') #disp('to the ACF in order to avoid this.') @@ -425,57 +425,57 @@ class CovData1D(WafoData): S[k] = 0. - ix = flatnonzero(k>2*I) - if ix.size>0: + ix = flatnonzero(k > 2 * I) + if ix.size > 0: ## # truncating all oscillating values above 2 times the peak ## # frequency to zero to ensure that ## # that high frequency noise is not added to ## # the simulated timeseries. ix0 = k[ix[0]] - S[ix0:-ix0] =0.0 + S[ix0:-ix0] = 0.0 trunc = 1e-5 maxS = S[I] - k = flatnonzero(S[I:-I]0: - S[k+I]=0. + k = flatnonzero(S[I:-I] < maxS * trunc) + if k.size > 0: + S[k + I] = 0. ## truncating small values to zero to ensure that ## that high frequency noise is not added to ## the simulated timeseries - cases1 = floor(cases/2) - cases2 = ceil(cases/2) + cases1 = floor(cases / 2) + cases2 = ceil(cases / 2) # Generate standard normal random numbers for the simulations #randn = np.random.randn - epsi = randn(nfft,cases2)+1j*randn(nfft,cases2) - Ssqr = sqrt(S/(nfft)) # #sqrt(S(wn)*dw ) - ephat = epsi*Ssqr #[:,np.newaxis] - y = fft(ephat,nfft,axis=0) - x[:, 1:cases+1] = hstack((y[2:ns+2, 0:cases2].real, y[2:ns+2, 0:cases1].imag)) + epsi = randn(nfft, cases2) + 1j * randn(nfft, cases2) + Ssqr = sqrt(S / (nfft)) # #sqrt(S(wn)*dw ) + ephat = epsi * Ssqr #[:,np.newaxis] + y = fft(ephat, nfft, axis=0) + x[:, 1:cases + 1] = hstack((y[2:ns + 2, 0:cases2].real, y[2:ns + 2, 0:cases1].imag)) - x[:, 0] = linspace(0,(ns-1)*dT,ns) ##(0:dT:(dT*(np-1)))' + x[:, 0] = linspace(0, (ns - 1) * dT, ns) ##(0:dT:(dT*(np-1)))' if derivative: - Ssqr = Ssqr*r_[0:(nfft/2+1), -(nfft/2-1):0]*2*pi/nfft/dT - ephat = epsi*Ssqr #[:,newaxis] - y = fft(ephat,nfft,axis=0) - xder[:, 1:(cases+1)] = hstack((y[2:ns+2, 0:cases2].imag -y[2:ns+2, 0:cases1].real)) - xder[:, 0] = x[:,0] + Ssqr = Ssqr * r_[0:(nfft / 2 + 1), -(nfft / 2 - 1):0] * 2 * pi / nfft / dT + ephat = epsi * Ssqr #[:,newaxis] + y = fft(ephat, nfft, axis=0) + xder[:, 1:(cases + 1)] = hstack((y[2:ns + 2, 0:cases2].imag - y[2:ns + 2, 0:cases1].real)) + xder[:, 0] = x[:, 0] if self.tr is not None: print(' Transforming data.') g = self.tr if derivative: for ix in range(cases): - tmp = g.gauss2dat(x[:,ix+1], xder[:,ix+1]) - x[:,ix+1] = tmp[0] - xder[:,ix+1] = tmp[1] + tmp = g.gauss2dat(x[:, ix + 1], xder[:, ix + 1]) + x[:, ix + 1] = tmp[0] + xder[:, ix + 1] = tmp[1] else: for ix in range(cases): - x[:, ix+1] = g.gauss2dat(x[:, ix+1]) + x[:, ix + 1] = g.gauss2dat(x[:, ix + 1]) if derivative: return x, xder @@ -574,7 +574,7 @@ class CovData1D(WafoData): warnings.warn(txt) return self.sim(ns=N, cases=cases), zeros(Ns), zeros(Ns) - indg = where(1-isnan(x))[0] #indices to the known observations + indg = where(1 - isnan(x))[0] #indices to the known observations #initializing variables mu1o = zeros(Ns, 1) @@ -588,121 +588,121 @@ class CovData1D(WafoData): if method.startswith('dec1'): # only correct for variables having the Markov property # but still seems to give a reasonable answer. Slow procedure. - Sigma = sptoeplitz(hstack((acf, zeros(N-n)))) + Sigma = sptoeplitz(hstack((acf, zeros(N - n)))) #Soo=Sigma(~inds,~inds); # covariance between known observations #S11=Sigma(inds,inds); # covariance between unknown observations #S1o=Sigma(inds,~inds);# covariance between known and unknown observations #tmp=S1o*pinv(full(Soo)); #tmp=S1o/Soo; # this is time consuming if Soo large - tmp = 2*Sigma[inds, indg]/(Sigma[indg, indg] + Sigma[indg, indg].T ) + tmp = 2 * Sigma[inds, indg] / (Sigma[indg, indg] + Sigma[indg, indg].T) if compute_sigma: #standard deviation of the expected surface #mu1o_std=sqrt(diag(S11-tmp*S1o')); - mu1o_std = sqrt(diag(Sigma[inds, inds]-tmp*Sigma[indg, inds])) + mu1o_std = sqrt(diag(Sigma[inds, inds] - tmp * Sigma[indg, inds])) #expected surface conditioned on the known observations from x - mu1o = tmp*x[indg] + mu1o = tmp * x[indg] #expected surface conditioned on the known observations from xs - mu1os = tmp*(xs[indg,:]) + mu1os = tmp * (xs[indg, :]) # sampled surface conditioned on the known observations - sample = mu1o + xs[inds,:] - mu1os + sample = mu1o + xs[inds, :] - mu1os elif method.startswith('dec2'): # only correct for variables having the Markov property # but still seems to give a reasonable answer # approximating the expected surfaces conditioned on # the known observations from x and xs by only using the closest points - Sigma = sptoeplitz(hstack((acf,zeros(n)))) - n2 = int(floor(n/2)) - idx = r_[0:2*n] + max(0,inds[0]-n2) # indices to the points used - tmpinds = zeros(N,dtype=bool) + Sigma = sptoeplitz(hstack((acf, zeros(n)))) + n2 = int(floor(n / 2)) + idx = r_[0:2 * n] + max(0, inds[0] - n2) # indices to the points used + tmpinds = zeros(N, dtype=bool) tmpinds[inds] = True # temporary storage of indices to missing points tinds = where(tmpinds[idx])[0] # indices to the points used - tindg = where(1-tmpinds[idx])[0] + tindg = where(1 - tmpinds[idx])[0] ns = len(tinds); # number of missing data in the interval nprev = 0; # number of previously simulated points - xsinds = xs[inds,:] - while ns>0: - tmp=2*Sigma[tinds, tindg]/(Sigma[tindg, tindg]+Sigma[tindg, tindg].T) + xsinds = xs[inds, :] + while ns > 0: + tmp = 2 * Sigma[tinds, tindg] / (Sigma[tindg, tindg] + Sigma[tindg, tindg].T) if compute_sigma: #standard deviation of the expected surface #mu1o_std=sqrt(diag(S11-tmp*S1o')); - ix = slice(nprev+1,nprev+ns+1) - mu1o_std[ix] = max(mu1o_std[ix], - sqrt(diag(Sigma[tinds, tinds]-tmp*Sigma[tindg,tinds]))) + ix = slice(nprev + 1, nprev + ns + 1) + mu1o_std[ix] = max(mu1o_std[ix], + sqrt(diag(Sigma[tinds, tinds] - tmp * Sigma[tindg, tinds]))) #end #expected surface conditioned on the closest known observations # from x and xs2 - mu1o[(nprev+1):(nprev+ns+1)] = tmp*x[idx[tindg]] - mu1os[(nprev+1):(nprev+ns+1),:] = tmp*xs[idx[tindg],:] + mu1o[(nprev + 1):(nprev + ns + 1)] = tmp * x[idx[tindg]] + mu1os[(nprev + 1):(nprev + ns + 1), :] = tmp * xs[idx[tindg], :] - if idx[-1]==N-1:# - ns =0 # no more points to simulate + if idx[-1] == N - 1:# + ns = 0 # no more points to simulate else: # updating by putting expected surface into x - x[idx[tinds]] = mu1o[(nprev+1):(nprev+ns+1)] - xs[idx[tinds]] = mu1os[(nprev+1):(nprev+ns+1)] + x[idx[tinds]] = mu1o[(nprev + 1):(nprev + ns + 1)] + xs[idx[tinds]] = mu1os[(nprev + 1):(nprev + ns + 1)] nw = sum(tmpinds[idx[-n2:]])# # data which we want to simulate once tmpinds[idx[:-n2]] = False # removing indices to data .. # which has been simulated - nprev = nprev+ns-nw # update # points simulated so far + nprev = nprev + ns - nw # update # points simulated so far - if (nw==0) and (nprev0: + while ns > 0: #make sure MATLAB uses a symmetric matrix solver - tmp = 2*Sigma[tinds,tindg]/(Sigma[tindg,tindg]+Sigma[tindg,tindg].T) - Sigma1o = Sigma[tinds,tinds] - tmp*Sigma[tindg,tinds] + tmp = 2 * Sigma[tinds, tindg] / (Sigma[tindg, tindg] + Sigma[tindg, tindg].T) + Sigma1o = Sigma[tinds, tinds] - tmp * Sigma[tindg, tinds] if compute_sigma: #standard deviation of the expected surface #mu1o_std=sqrt(diag(S11-tmp*S1o')); - mu1o_std[(nprev+1):(nprev+ns+1)] = max( mu1o_std[(nprev+1):(nprev+ns)] , + mu1o_std[(nprev + 1):(nprev + ns + 1)] = max(mu1o_std[(nprev + 1):(nprev + ns)] , sqrt(diag(Sigma1o))) #end #expected surface conditioned on the closest known observations from x - mu1o[(nprev+1):(nprev+ns+1)] = tmp*x2[idx[tindg]] + mu1o[(nprev + 1):(nprev + ns + 1)] = tmp * x2[idx[tindg]] #sample conditioned on the known observations from x - sample[(nprev+1):(nprev+ns+1),:] = rndnormnd(tmp*x[idx[tindg]],Sigma1o, cases) - if idx[-1] == N-1: + sample[(nprev + 1):(nprev + ns + 1), :] = rndnormnd(tmp * x[idx[tindg]], Sigma1o, cases) + if idx[-1] == N - 1: ns = 0 # no more points to simulate else: # updating - x2[idx[tinds]] = mu1o[(nprev+1):(nprev+ns+1)] #expected surface - x[idx[tinds]] = sample[(nprev+1):(nprev+ns+1)]#sampled surface - nw = sum(tmpinds[idx[-n2::]]==True)# # data we want to simulate once more + x2[idx[tinds]] = mu1o[(nprev + 1):(nprev + ns + 1)] #expected surface + x[idx[tinds]] = sample[(nprev + 1):(nprev + ns + 1)]#sampled surface + nw = sum(tmpinds[idx[-n2::]] == True)# # data we want to simulate once more tmpinds[idx[:-n2]] = False # removing indices to data .. # which has been simulated - nprev = nprev+ns-nw # update # points simulated so far + nprev = nprev + ns - nw # update # points simulated so far - if (nw==0) and (nprev0.3*n: + if len(k) > 0.3 * n: return toeplitz(x) else: spdiags = sparse.dia_matrix - data = x[k].reshape(-1,1).repeat(n,axis=-1) + data = x[k].reshape(-1, 1).repeat(n, axis= -1) offsets = k - y = spdiags((data, offsets), shape=(n,n)) - if k[0]==0: + y = spdiags((data, offsets), shape=(n, n)) + if k[0] == 0: offsets = k[1::] - data = data[1::,:] - return y + spdiags((data, -offsets), shape=(n,n)) + data = data[1::, :] + return y + spdiags((data, -offsets), shape=(n, n)) def test_covdata(): import wafo.data @@ -823,4 +823,4 @@ if __name__ == '__main__': import doctest doctest.testmod() else: - main() \ No newline at end of file + main() diff --git a/pywafo/src/wafo/data_structures.py b/pywafo/src/wafo/data_structures.py deleted file mode 100644 index 23d59dc..0000000 --- a/pywafo/src/wafo/data_structures.py +++ /dev/null @@ -1,2659 +0,0 @@ -""" -Class objects for data, spectrum, covariance function and time series in WAFO - - To represent data, spectra, covariance functions and time series - in WAFO, the classes WafoData, SpecData1D, CovData1D and TimeSeries is used. - Here follows a list of their attributes: - - WafoData - ---------- - data : array_like - args : vector for 1D, list of vectors for 2D, 3D, ... - date : Date and time of creation or change. - labels : AxisLabels - children : list of WafoData objects - - SpecData1D - ------------ - data : One sided Spectrum values - args : freguency values of freqtype - type :String: 'freq', 'dir', 'k2d', k1d', 'encdir' or 'enc'. - freqtype :'w' OR 'f' OR 'k' Frequency/wave number lag, length nf. - tr : Transformation function (default [] (none)). - h : Water depth (default inf). - norm : Normalization flag, Logical 1 if S is normalized, 0 if not - date :Date and time of creation or change. - v : Ship speed, if .type = 'enc' or 'encdir'. - phi angle of rotation of the coordinate system - (counter-clocwise) e.g. azymuth of a ship. - - - CovData1D - ---------- - data : Covariance function values. Size [ny nx nt], all singleton dim. removed. - args : Lag of first space dimension, length nx. - .h Water depth. - .tr Transformation function. - .type 'enc', 'rot' or 'none'. - .v Ship speed, if .type='enc' - .phi Rotation of coordinate system, e.g. direction of ship - .norm Normalization flag, Logical 1 if autocorrelation, 0 if covariance. - .Rx ... .Rtttt Obvious derivatives of .R. - .note Memorandum string. - .date Date and time of creation or change. - - - - -""" -#------------------------------------------------------------------------------- -# Name: module1 -# Purpose: -# -# Author: pab -# -# Created: 16.09.2008 -# Copyright: (c) pab 2008 -# Licence: -#------------------------------------------------------------------------------- -#!/usr/bin/env python - - -from __future__ import division -import warnings -import numpy -import numpy as np -from time import gmtime, strftime -from wafo.misc import nextpow2, discretize, findextrema, findrfc, sub_dict_select -#import wafo.spectrum.dispersion_relation as disp_rel -#from scipy import fft -from scipy.integrate import simps, trapz -from pylab import stineman_interp -import scipy.interpolate as interpolate - -import warnings -try: - import diffsumfunq -except: - pass - -from plotbackend import plotbackend - -__all__ = ['SpecData1D','SpecData2D','WafoData', 'AxisLabels','CovData1D', - 'TimeSeries','sensortypeid','sensortype'] -def empty_copy(obj): - class Empty(obj.__class__): - def __init__(self): - pass - newcopy = Empty() - newcopy.__class__ = obj.__class__ - return newcopy - -def _set_seed(iseed): - if iseed != None: - try: - np.random.set_state(iseed) - except: - np.random.seed(iseed) - -class WafoData(object): - '''Container class for data objects in WAFO - - Member variables - ---------------- - data : array_like - args : vector for 1D, list of vectors for 2D, 3D, ... - labels : AxisLabels - children : list of WafoData objects - - Member methods - -------------- - plot : - copy : - - - Example - ------- - >>> import numpy as np - >>> x = np.arange(-2,2,0.2) - - # Plot 2 objects in one call - >>> d2 = WafoData(np.sin(x),x,xlab='x',ylab='sin',title='sinus') - >>> d2.plot() - - % Plot with confidence interval - d3 = wdata(sin(x),x); - d3 = set(d3,'dataCI',[sin(x(:))*0.9 sin(x(:))*1.2]); - plot(d3) - - See also - -------- - wdata/plot, - specdata, - covdata - ''' - def __init__(self,data=None,args=None,**kwds): - self.data = data - self.args = args - self.date = strftime("%a, %d %b %Y %H:%M:%S", gmtime()) - self.plotter = None - self.children = None - self.labels = AxisLabels(**kwds) - self.setplotter() - - def plot(self,*args,**kwds): - if self.children!=None: - plotbackend.hold('on') - for child in self.children: - tmp = child.plot(*args,**kwds) - - tmp2 = self.plotter.plot(self,*args,**kwds) - return tmp2 - - def copy(self): - newcopy = empty_copy(self) - newcopy.__dict__.update(self.__dict__) - return newcopy - - - def setplotter(self): - ''' - Set plotter based on the data type data_1d, data_2d, data_3d or data_nd - ''' - - if isinstance(self.args,(list,tuple)): # Multidimensional data - ndim = len(self.args) - if ndim<2: - warnings.warn('Unable to determine plotter-type, because len(self.args)<2.') - print('If the data is 1D, then self.args should be a vector!') - print('If the data is 2D, then length(self.args) should be 2.') - print('If the data is 3D, then length(self.args) should be 3.') - print('Unless you fix this, the plot methods will not work!') - elif ndim==2: - self.plotter = Plotter_2d() - else: - warnings.warn('Plotter method not implemented for ndim>2') - - else: #One dimensional data - self.plotter = Plotter_1d() - - -class AxisLabels: - def __init__(self,title='',xlab='',ylab='',zlab='',**kwds): - self.title = title - self.xlab = xlab - self.ylab = ylab - self.zlab = zlab - def copy(self): - lbkwds = self.labels.__dict__.copy() - labels = AxisLabels(**lbkwds) - return labels - - def labelfig(self): - try: - plotbackend.title(self.title) - plotbackend.xlabel(self.xlab) - plotbackend.ylabel(self.ylab) - plotbackend.zlabel(self.zlab) - except: - pass - -class Plotter_1d(object): - """ - class comment - - bar - barh - loglog - semilogx - semilogy - plot - stem - scatter - - """ - - def __init__(self,plotmethod='plot'): - self.plotfun = None - self.plotbackend = plotbackend - try: - self.plotfun = plotbackend.__dict__[plotmethod] - except: - pass - def show(self): - plotbackend.show() - - def plot(self,wdata,*args,**kwds): - if isinstance(wdata.args,(list,tuple)): - args1 = tuple((wdata.args))+(wdata.data,)+args - else: - args1 = tuple((wdata.args,))+(wdata.data,)+args - self.plotfun(*args1,**kwds) - wdata.labels.labelfig() - -class Plotter_2d(Plotter_1d): - """ - class comment - - contour - mesh - surf - - - """ - - def __init__(self,plotmethod='contour'): - super(Plotter_2d,self).__init__(plotmethod) - #self.plotfun = plotbackend.__dict__[plotmethod] - -class TrGauss(WafoData): - def __init__(self,*args,**kwds): - super(TrGauss, self).__init__(*args,**kwds) - -class CycleCount(object): - def __init__(self,M,m): - self.M = M - self.m = m - - def amplitudes(self): - return (self.M-self.m)/2. - -class SpecData1D(WafoData): - """ Container class for 1D spectrum data objects in WAFO - - Member variables - ---------------- - data : array_like - args : vector for 1D, list of vectors for 2D, 3D, ... - - type : string - spectrum type, one of 'freq', 'k1d', 'enc' (default 'freq') - freqtype : letter - frequency type, one of: 'f', 'w' or 'k' (default 'w') - - - Examples - -------- - >>> import numpy as np - >>> import wafo.spectrum.models as sm - >>> Sj = sm.Jonswap(Hm0=3) - >>> w = np.linspace(0,4,256) - >>> S1 = Sj.toSpecData(w) #Make spectrum object from numerical values - >>> S = SpecData1D(Sj(w),w) # Alternatively do it manually - - See also - -------- - WafoData - CovData - """ - - def __init__(self,*args,**kwds): - super(SpecData1D, self).__init__(*args,**kwds) - - self.name='WAFO Spectrum Object' - self.type='freq' - self.freqtype='w' - self.angletype='' - self.h=np.inf - self.tr=None - self.phi=0. - self.v=0. - self.norm=0 - somekeys = ['angletype', 'phi', 'name', 'h', 'tr', 'freqtype', 'v', 'type', 'norm'] - - self.__dict__.update(sub_dict_select(kwds,somekeys)) - - self.setlabels() -## def copy(self): -## newcopy = empty_copy(self) -## newcopy.update(self.__dict__) -## #kwds = self.__dict__.copy() -## #wdata = SpecData1D(**kwds) -## #wdata.labels = sel.labels.copy() -## return newcopy - def toacf_matrix(self,nr=0,Nt=None,dt=None): - ''' Computes covariance function and its derivatives, alternative version - - Parameters - ---------- - Nr = number of derivatives in output, Nr<=4 (default 0) - Nt = number in time grid, i.e., number of time-lags. - (default rate*(n-1)) where rate = round(1/(2*f(end)*dt)) or - rate = round(pi/(w(n)*dt)) depending on S. - dt = time spacing for R - - Returns - ------- - R = [R0, R1,...Rnr] matrix with autocovariance and its - derivatives, i.e., Ri (i=1:nr) are column vectors with - the 1'st to nr'th derivatives of R0. size Nt+1 x Nr+1 - - NB! This routine requires that the spectrum grid is equidistant - starting from zero frequency. - Example: - -------- - S = jonswap; - dt = 0.1; - R = spec2cov2(S,3,256,dt); - See also - -------- - spec2cov, specinterp, datastructures - ''' - - ftype = self.freqtype #; %options are 'f' and 'w' and 'k' - freq = self.args - n = length(freq) - dTold = self.sampling_period() - if dt is None: - dt = dTold - rate = 1 - else: - rate = np.maximum(np.round(dTold*1./dt),1.) - - - if Nt is None: - Nt = rate*(n-1) - else: #%check if Nt is ok - Nt = np.minimum(Nt,rate*(n-1)); - - - checkdt = 1.2*min(np.diff(freq))/2./np.pi; - if ftype in 'k': - lagtype = 'x' - else: - lagtype = 't' - if ftype in 'f': - checkdt = checkdt*2*np.pi - msg1 = 'The step dt = %g in computation of the density is too small.' % dt - msg2 = 'The step dt = %g step is small, may cause numerical inaccuracies.' % dt - - if (checkdt < 2.**-16/dt): - np.disp(msg1) - np.disp('The computed covariance (by FFT(2^K)) may differ from the theoretical.') - np.disp('Solution:') - raise ValueError('use larger dt or sparser grid for spectrum.') - - - #% Calculating covariances - #%~~~~~~~~~~~~~~~~~~~~~~~~ - S2 = self.copy() - S2.resample(dt) - - R2 = S2.toacf(nr,Nt,rate=1) - R = np.zeros((Nt+1,nr+1)) - R[:,0] = R2.data[0:Nt+1] - fieldname = 'R' + lagtype*nr - for ix in range(1,nr+1): - fn = fieldname[1:ix]; - R[:,0:Nt+1] = gettattr(R2,fn)[0:Nt+1]; - - - EPS0 = 0.0001 - cc = R[0,0]-R[1,0]*(R[1,0]/R[0,0]); - if Nt+1>=5: - #%cc1=R(1,1)-R(3,1)*(R(3,1)/R(1,1))+R(3,2)*(R(3,2)/R(1,3)); - #%cc3=R(1,1)-R(5,1)*(R(5,1)/R(1,1))+R(5,2)*(R(5,2)/R(1,3)); - - cc2 = R(0,0)-R(4,0)*(R(4,0)/R(0,0)); - if (cc2R(x,t); Nt and Ny set =>R(y,t) - 4) Any type, Nt, Nx and Ny set => R(x,y,t) - 5) Any type, Nt not set, Nx and/or Ny set => Nt set to default, goto 3) or 4) - - NB! This routine requires that the spectrum grid is equidistant - starting from zero frequency. - NB! If you are using a model spectrum, S, with sharp edges - to calculate covariances then you should probably round off the sharp - edges like this: - - Example: - >>> import wafo.spectrum.models as sm - >>> Sj = sm.Jonswap() - >>> S = Sj.toSpecData() - >>> S.data[0:40] = 0.0 - >>> S.data[100:-1] = 0.0 - >>> Nt = len(S.data)-1; - >>> R = S.toacf(nr=0,Nt=Nt) - - R = spec2cov(S,0,Nt); - win = parzen(2*Nt+1); - R.R = R.R.*win(Nt+1:end); - S1 = cov2spec(R); - R2 = spec2cov(S1); - figure(1) - plotspec(S),hold on, plotspec(S1,'r') - figure(2) - covplot(R), hold on, covplot(R2,[],[],'r') - figure(3) - semilogy(abs(R2.R-R.R)), hold on, - semilogy(abs(S1.S-S.S)+1e-7,'r') - - See also cov2spec, datastructures - ''' - - freq = self.args - n = freq.size - - if freq[0]>0: - raise ValueError('Spectrum does not start at zero frequency/wave number.\n Correct it with resample, for example.') - dw = np.abs(np.diff(freq,n=2,axis=0)) - if np.any(dw>1.0e-8): - raise ValueError('Not equidistant frequencies/wave numbers in spectrum.\n Correct it with resample, for example.') - - - if rate is None: - rate=1 #; %interpolation rate - elif rate>16: - rate=16 - else: # make sure rate is a power of 2 - rate=2**nextpow2(rate) - - if Nt is None: - Nt = rate*(n-1) - else: #check if Nt is ok - Nt = np.minimum(Nt,rate*(n-1)) - - S = self.copy() - - - if self.freqtype in 'k': - lagtype = 'x' - else: - lagtype = 't' - - dT = S.sampling_period() - #normalize spec so that sum(specn)/(n-1)=R(0)=var(X) - specn = S.data*freq[-1] - if S.freqtype in 'f': - w = freq*2*np.pi - else: - w = freq - - nfft = rate*2**nextpow2(2*n-2); - - Rper = np.r_[specn, np.zeros(nfft-(2*n)+2) , np.conj(specn[n-1:0:-1])] #; % periodogram - t = np.r_[0:Nt+1]*dT*(2*n-2)/nfft - - fft = np.fft.fft - - r = fft(Rper,nfft).real/(2*n-2) - R = CovData1D(r[0:Nt+1],t,lagtype=lagtype) - R.tr = S.tr; - R.h = S.h; - R.norm = S.norm; - - if nr>0: - w = np.r_[w , np.zeros(nfft-2*n+2) ,-w[n-1:0:-1] ] - fieldname = 'R' + lagtype*nr ; - for ix in range(1,nr+1): - Rper = (-1j*w*Rper) - r = fft(Rper,nfft).real/(2*n-2) - setattr(R,fieldname[0:ix+1],r[0:Nt+1]) - return R - - - - - def sim(self,ns=None,cases=1,dt=None,iseed=None,method='random',derivative=False): - ''' Simulates a Gaussian process and its derivative from spectrum - - Parameters - ---------- - ns : scalar - number of simulated points. (default length(S)-1=n-1). - If ns>n-1 it is assummed that R(k)=0 for all k>n-1 - cases : scalar - number of replicates (default=1) - dt : scalar - step in grid (default dt is defined by the Nyquist freq) - iseed : int or state - starting state/seed number for the random number generator - (default none is set) - method : string - if 'exact' : simulation using cov2sdat - if 'random' : random phase and amplitude simulation (default) - derivative : bool - if true : return derivative of simulated signal as well - otherwise - - - Returns - ------- - xs = a cases+1 column matrix ( t,X1(t) X2(t) ...). - xsder = a cases+1 column matrix ( t,X1'(t) X2'(t) ...). - - Details - ------- - Performs a fast and exact simulation of stationary zero mean - Gaussian process through circulant embedding of the covariance matrix - or by summation of sinus functions with random amplitudes and random - phase angle. - - If the spectrum has a non-empty field .tr, then the transformation is - applied to the simulated data, the result is a simulation of a transformed - Gaussian process. - - Note: The method 'exact' simulation may give high frequency ripple when - used with a small dt. In this case the method 'random' works better. - - Example: - >>> import wafo.spectrum.models as sm - >>> Sj = sm.Jonswap();S = Sj.toSpecData(); - >>> np =100; dt = .2; - >>> x1 = S.sim(np,dt=dt); - - waveplot(x1,'r',x2,'g',1,1) - - See also - -------- - cov2sdat, gaus2dat - - Reference - ----------- - C.R Dietrich and G. N. Newsam (1997) - "Fast and exact simulation of stationary - Gaussian process through circulant embedding - of the Covariance matrix" - SIAM J. SCI. COMPT. Vol 18, No 4, pp. 1088-1107 - - Hudspeth, R.T. and Borgman, L.E. (1979) - "Efficient FFT simulation of Digital Time sequences" - Journal of the Engineering Mechanics Division, ASCE, Vol. 105, No. EM2, - - ''' - - fft = np.fft.fft - - S = self.copy() - if dt is not None: - S.resample(dt) - - - ftype = S.freqtype - freq = S.args - - dT = S.sampling_period() - Nt = freq.size - - if ns is None: - ns=Nt-1 - - if method in 'exact': - - #nr=0,Nt=None,dt=None - R = S.toacf(nr=0); - T = Nt*dT - ix = np.flatnonzero(R.args>T); - - # Trick to avoid adding high frequency noise to the spectrum - if ix.size>0: - R.data[ix[0]::]=0.0 - - return R.sim(ns=ns,cases=cases,iseed=iseed,derivative=derivative) - - _set_seed(iseed) - - ns = ns+np.mod(ns,2) # make sure it is even - - fi = freq[1:-1] - Si = S.data[1:-1] - if ftype in ('w','k'): - fact = 2.*np.pi - Si = Si*fact; - fi = fi/fact; - - zeros = np.zeros - - x = zeros((ns,cases+1)); - - df = 1/(ns*dT) - - - # interpolate for freq. [1:(N/2)-1]*df and create 2-sided, uncentered spectra - # ---------------------------------------------------------------------------- - f = np.arange(1,ns/2.)*df - - Fs = np.hstack((0., fi, df*ns/2.)) - Su = np.hstack((0., np.abs(Si)/2., 0.)); - - - Si = np.interp(f,Fs,Su) - Su=np.hstack((0., Si,0, Si[(ns/2)-2::-1])) - del(Si, Fs) - - # Generate standard normal random numbers for the simulations - # ----------------------------------------------------------- - randn = np.random.randn - Zr = randn((ns/2)+1,cases) - Zi = np.vstack((zeros((1,cases)), randn((ns/2)-1,cases), zeros((1,cases)))); - - A = zeros((ns,cases),dtype=complex) - A[0:(ns/2+1),:] = Zr - 1j*Zi; - del(Zr, Zi) - A[(ns/2+1):ns,:] = A[ns/2-1:0:-1,:].conj(); - A[0,:] = A[0,:]*np.sqrt(2.); - A[(ns/2),:] = A[(ns/2),:]*np.sqrt(2.); - - - # Make simulated time series - # -------------------------- - - T = (ns-1)*dT - Ssqr = np.sqrt(Su*df/2.) - - # stochastic amplitude - A = A*Ssqr[:,np.newaxis]; - - - # Deterministic amplitude - #A = sqrt[1]*Ssqr(:,ones(1,cases)).*exp(sqrt(-1)*atan2(imag(A),real(A))); - del( Su, Ssqr) - - - x[:,1::] = fft(A,axis=0).real - x[:,0] = np.linspace(0,T,ns) #'; %(0:dT:(np-1)*dT).'; - - - if derivative: - xder=np.zeros(ns,cases+1) - w = 2.*np.pi*np.hstack((0, f, 0.,-f[-1::-1])) - A = -1j*A*w[:,newaxis] - xder[:,1:(cases+1)] = fft(A,axis=0).real; - xder[:,0] = x[:,0] - - - - if S.tr is not None: - np.disp(' Transforming data.') - g=S.tr; - G=np.fliplr(g); #% the invers of g - if derivative: - for ix in range(cases): - tmp=tranproc(np.hstack((x[:,ix+1], xder[:,ix+1])),G); - x[:,ix+1]=tmp[:,0]; - xder[:,ix+1]=tmp[:,1]; - - else: - for ix in range(cases): - x[:,ix+1]=tranproc(x[:,ix+1],G); - - - - if derivative: - return x,xder - else: - return x - -# function [x2,x,svec,dvec,A]=spec2nlsdat(S,np,dt,iseed,method,truncationLimit) - def sim_nl(self,ns=None,cases=1,dt=None,iseed=None,method='random', - fnlimit=1.4142,reltol=1e-3,g=9.81): - """ Simulates a Randomized 2nd order non-linear wave X(t) - - Parameters - ---------- - ns : scalar - number of simulated points. (default length(S)-1=n-1). - If ns>n-1 it is assummed that R(k)=0 for all k>n-1 - cases : scalar - number of replicates (default=1) - dt : scalar - step in grid (default dt is defined by the Nyquist freq) - iseed : int or state - starting state/seed number for the random number generator - (default none is set) - method : string - 'apStochastic' : Random amplitude and phase (default) - 'aDeterministic' : Deterministic amplitude and random phase - 'apDeterministic' : Deterministic amplitude and phase - fnLimit : scalar - normalized upper frequency limit of spectrum for 2'nd order - components. The frequency is normalized with - sqrt(gravity*tanh(kbar*waterDepth)/Amax)/(2*pi) - (default sqrt(2), i.e., Convergence criterion [1]_). - Other possible values are: - sqrt(1/2) : No bump in trough criterion - sqrt(pi/7) : Wave steepness criterion - reltol : scalar - relative tolerance defining where to truncate spectrum for the - for sum and difference frequency effects - - - Returns - ------- - xs2 = a cases+1 column matrix ( t,X1(t) X2(t) ...). - xs1 = a cases+1 column matrix ( t,X1'(t) X2'(t) ...). - - Details - ------- - Performs a Fast simulation of Randomized 2nd order non-linear - waves by summation of sinus functions with random amplitudes and - phase angles. The extent to which the simulated result are applicable - to real seastates are dependent on the validity of the assumptions: - - 1. Seastate is unidirectional - 2. Surface elevation is adequately represented by 2nd order random - wave theory - 3. The first order component of the surface elevation is a Gaussian - random process. - - If the spectrum does not decay rapidly enough towards zero, the - contribution from the 2nd order wave components at the upper tail can - be very large and unphysical. To ensure convergence of the perturbation - series, the upper tail of the spectrum is truncated at FNLIMIT in the - calculation of the 2nd order wave components, i.e., in the calculation - of sum and difference frequency effects. This may also be combined with - the elimination of second order effects from the spectrum, i.e., extract - the linear components from the spectrum. One way to do this is to use - SPEC2LINSPEC. - - Example - -------- - np =100; dt = .2; - [x1, x2] = spec2nlsdat(jonswap,np,dt); - waveplot(x1,'r',x2,'g',1,1) - - See also - -------- - spec2linspec, spec2sdat, cov2sdat - - References - ---------- - .. [1] Nestegaard, A and Stokka T (1995) - A Third Order Random Wave model. - In proc.ISOPE conf., Vol III, pp 136-142. - - .. [2] R. S Langley (1987) - A statistical analysis of non-linear random waves. - Ocean Engng, Vol 14, pp 389-407 - - .. [3] Marthinsen, T. and Winterstein, S.R (1992) - 'On the skewness of random surface waves' - In proc. ISOPE Conf., San Francisco, 14-19 june. - """ - - - # TODO % Check the methods: 'apdeterministic' and 'adeterministic' - - from wafo.spectrum import dispersion_relation as sm - - Hm0 = self.characteristic('Hm0')[0] - Tm02 = self.characteristic('Tm02')[0] - #Hm0 = spec2char(S,'Hm0'); - #Tm02 = spec2char(S,'Tm02'); - - _set_seed(iseed) - fft = np.fft.fft - - S = self.copy() - if dt is not None: - S.resample(dt) - - - ftype = S.freqtype - freq = S.args - - dT = S.sampling_period() - Nt = freq.size - - if ns is None: - ns = Nt-1 - - ns = ns+np.mod(ns,2) # make sure it is even - - fi = freq[1:-1] - Si = S.data[1:-1] - if ftype in ('w','k'): - fact = 2.*np.pi - Si = Si*fact; - fi = fi/fact; - - Smax = max(Si) - waterDepth = min(abs(S.h),10.**30); - - zeros = np.zeros - - x = zeros((ns,cases+1)); - - df = 1/(ns*dT) - - # interpolate for freq. [1:(N/2)-1]*df and create 2-sided, uncentered spectra - # ---------------------------------------------------------------------------- - f = np.arange(1,ns/2.)*df - Fs = np.hstack((0., fi, df*ns/2.)) - w = 2.*np.pi*Fs - kw = sm.w2k(w ,0.,waterDepth,g)[0] - Su = np.hstack((0., np.abs(Si)/2., 0.)); - - - - Si = np.interp(f,Fs,Su) - nmin = (Si>Smax*reltol).argmax() - nmax = np.flatnonzero(Si>0).max() - Su = np.hstack((0., Si,0, Si[(ns/2)-2::-1])) - del(Si, Fs) - - # Generate standard normal random numbers for the simulations - # ----------------------------------------------------------- - randn = np.random.randn - Zr = randn((ns/2)+1,cases) - Zi = np.vstack((zeros((1,cases)), randn((ns/2)-1,cases), zeros((1,cases)))); - - A = zeros((ns,cases),dtype=complex) - A[0:(ns/2+1),:] = Zr - 1j*Zi; - del(Zr, Zi) - A[(ns/2+1):ns,:] = A[ns/2-1:0:-1,:].conj(); - A[0,:] = A[0,:]*np.sqrt(2.); - A[(ns/2),:] = A[(ns/2),:]*np.sqrt(2.); - - - # Make simulated time series - # -------------------------- - - T = (ns-1)*dT - Ssqr = np.sqrt(Su*df/2.) - - - if method.startswith('apd') : # apdeterministic - # Deterministic amplitude and phase - A[1:(ns/2),:] = A[1,0]; - A[(ns/2+1):ns,:] = A[1,0].conj(); - A = sqrt(2)*Ssqr[:,np.newaxis]*exp(1J*atan2(A.imag,A.real)) - elif method.startswith('ade'): # adeterministic - # Deterministic amplitude and random phase - A = sqrt(2)*Ssqr[:,np.newaxis]*exp(1J*atan2(A.imag,A.real)); - else: - # stochastic amplitude - A = A*Ssqr[:,np.newaxis]; - # Deterministic amplitude - #A = sqrt(2)*Ssqr(:,ones(1,cases)).*exp(sqrt(-1)*atan2(imag(A),real(A))); - del( Su, Ssqr) - - - x[:,1::] = fft(A,axis=0).real - x[:,0] = np.linspace(0,T,ns) #'; %(0:dT:(np-1)*dT).'; - - - - x2 = x.copy() - - # If the spectrum does not decay rapidly enough towards zero, the - # contribution from the wave components at the upper tail can be very - # large and unphysical. - # To ensure convergence of the perturbation series, the upper tail of the - # spectrum is truncated in the calculation of sum and difference - # frequency effects. - # Find the critical wave frequency to ensure convergence. - - sqrt = np.sqrt - log = np.log - pi = np.pi - tanh = np.tanh - numWaves = 1000. # Typical number of waves in 3 hour seastate - kbar = sm.w2k(2.*np.pi/Tm02,0.,waterDepth)[0]; - Amax = sqrt(2*log(numWaves))*Hm0/4; #% Expected maximum amplitude for 1000 waves seastate - - fLimitUp = fnlimit*sqrt(g*tanh(kbar*waterDepth)/Amax)/(2*pi); - fLimitLo = sqrt(g*tanh(kbar*waterDepth)*Amax/waterDepth)/(2*pi*waterDepth); - - nmax = min(np.flatnonzero(f<=fLimitUp).max(),nmax)+1; - nmin = max(np.flatnonzero(fLimitLo<=f).min(),nmin)+1; - - #if isempty(nmax),nmax = np/2;end - #if isempty(nmin),nmin = 2;end % Must always be greater than 1 - fLimitUp = df*nmax; - fLimitLo = df*nmin; - - print('2nd order frequency Limits = %g,%g'% (fLimitLo, fLimitUp)) - - - -## if nargout>3, -## %compute the sum and frequency effects separately -## [svec, dvec] = disufq((A.'),w,kw,min(h,10^30),g,nmin,nmax); -## svec = svec.'; -## dvec = dvec.'; -## -## x2s = fft(svec); % 2'nd order sum frequency component -## x2d = fft(dvec); % 2'nd order difference frequency component -## -## % 1'st order + 2'nd order component. -## x2(:,2:end) =x(:,2:end)+ real(x2s(1:np,:))+real(x2d(1:np,:)); -## else - rvec,ivec = diffsumfunq.disufq(A.real,A.imag,w,kw,waterDepth,g,nmin,nmax) - - svec = rvec + 1J*ivec - x2o = fft(svec) # 2'nd order component - - - # 1'st order + 2'nd order component. - x2[:,1::] = x[:,1::]+ x2o[0:ns,:].real(); - - return x2, x - - - - - - - def moment(self,nr=2,even=True,j=0): - ''' Calculates spectral moments from spectrum - - Parameters - ---------- - nr : int - order of moments (recomended maximum 4) - even : bool - False for all moments, - True for only even orders - j : int - 0 or 1 - - Returns - ------- - m : list of moments - mtext : list of strings describing the elements of m, see below - - Details - ------- - Calculates spectral moments of up to order NR by use of - Simpson-integration. - - / / - mj_t^i = | w^i S(w)^(j+1) dw, or mj_x^i = | k^i S(k)^(j+1) dk - / / - - where k=w^2/gravity, i=0,1,...,NR - - The strings in output mtext have the same position in the list - as the corresponding numerical value has in output m - Notation in mtext: 'm0' is the variance, - 'm0x' is the first-order moment in x, - 'm0xx' is the second-order moment in x, - 'm0t' is the first-order moment in t, - etc. - For the calculation of moments see Baxevani et al. - - Example: - >>> import numpy as np - >>> import wafo.spectrum.models as sm - >>> Sj = sm.Jonswap(Hm0=3) - >>> w = np.linspace(0,4,256) - >>> S = SpecData1D(Sj(w),w) #Make spectrum object from numerical values - >>> S.moment() - ([0.56220770033914191, 0.35433180985851975], ['m0', 'm0tt']) - - References - ---------- - Baxevani A. et al. (2001) - Velocities for Random Surfaces - ''' - pi= np.pi - one_dim_spectra = ['freq','enc','k1d'] - if self.type not in one_dim_spectra: - raise ValueError('Unknown spectrum type!') - - f = np.ravel(self.args) - S = np.ravel(self.data) - if self.freqtype in ['f','w']: - vari = 't' - if self.freqtype=='f': - f = 2.*pi*f - S = S/(2.*pi) - else: - vari = 'x' - S1=np.abs(S)**(j+1.) - m = [simps(S1,x=f)] - mtxt = 'm%d' % j - mtext = [mtxt] - step = np.mod(even,2)+1 - df = f**step - for i in range(step,nr+1,step): - S1 = S1*df - m.append(simps(S1,x=f)) - mtext.append(mtxt+vari*i) - return m, mtext - - def sampling_period(self): - ''' Returns sampling interval from Nyquist frequency of spectrum - - Returns - --------- - dT : scalar - sampling interval, unit: - [m] if wave number spectrum, - [s] otherwise - - Let wm be maximum frequency/wave number in spectrum, - then dT=pi/wm if angular frequency, dT=1/(2*wm) if natural frequency (Hz) - - Example - ------- - S = jonswap; - dt = spec2dt(S) - - See also - ''' - - if self.freqtype in 'f': - wmdt = 0.5 # Nyquist to sampling interval factor - else: # ftype == w og ftype == k - wmdt = np.pi; - - wm = self.args[-1] #Nyquist frequency - dt = wmdt/wm; #sampling interval = 1/Fs - return dt - - def resample(self,dt=None,Nmin=0,Nmax=2**13+1,method='stineman'): - ''' Interpolate and zero-padd spectrum to change Nyquist freq. - - Parameters - ---------- - dt : scalar - wanted sampling interval (default as given by S, see spec2dt) - unit: [s] if frequency-spectrum, [m] if wave number spectrum - Nmin : scalar - minimum number of frequencies. - Nmax : scalar - minimum number of frequencies - method : string - interpolation method (options are 'linear', 'cubic' or 'stineman') - - To be used before simulation (e.g. spec2sdat) or evaluation of covariance - function (spec2cov) to directly get wanted sampling interval. - The input spectrum is interpolated and padded with zeros to reach - the right max-frequency, w(end)=pi/dt, f(end)=1/(2*dt), or k(end)=pi/dt. - The objective is that output frequency grid should be at least as dense - as the input grid, have equidistant spacing and length equal to - 2^k+1 (>=Nmin). If the max frequency is changed, the number of points - in the spectrum is maximized to 2^13+1. - - Note: Also zero-padding down to zero freq, if S does not start there. - If empty input dt, this is the only effect. - - See also - -------- - spec2cov, spec2sdat, covinterp, spec2dt - ''' - - - ftype = self.freqtype - w = self.args.ravel() - n = w.size - - #%doInterpolate = 0; - - if ftype=='f': - Cnf2dt = 0.5; # Nyquist to sampling interval factor - else: #% ftype == w og ftype == k - Cnf2dt = np.pi - - wnOld = w[-1] # Old Nyquist frequency - dTold = Cnf2dt/wnOld; # sampling interval=1/Fs - - - if dt is None: - dt=dTold - - - - # Find how many points that is needed - nfft = 2**nextpow2(max(n-1,Nmin-1)) - dttest = dTold*(n-1)/nfft - - while (dttest>dt) and (nfft0 or w[1]>0 or (nfft!=n) or dt!=dTold or any(abs(np.diff(w,axis=0))>1.0e-8); - - if doInterpolate>0: - S1 = self.data - - dw = min(np.diff(w)); - - if dWn>0: - #% add a zero just above old max-freq, and a zero at new max-freq - #% to get correct interpolation there - Nz = 1 + (dWn>dw) #; % Number of zeros to add - if Nz==2: - w = np.hstack((w, wnOld+dw, wnNew)) - else: - w = np.hstack((w, wnNew)) - - S1 = np.hstack((S1, np.zeros(Nz))) - - if w[0]>0: - #% add a zero at freq 0, and, if there is space, a zero just below min-freq - Nz = 1 + (w[0]>dw); #% Number of zeros to add - if Nz==2: - w=np.hstack((0, w[0]-dw, w)) - else: - w=np.hstack((0, w)) - - S1 = np.hstack((zeros(Nz), S1)) - - - #% Do a final check on spacing in order to check that the gridding is - #% sufficiently dense: - #np1 = S1.size - dwMin = np.finfo(float).max - #%wnc = min(wnNew,wnOld-1e-5); - wnc = wnNew; - specfun = lambda xi : stineman_interp(xi,w,S1) - - - x,y = discretize(specfun,0,wnc) - dwMin = np.minimum(min(np.diff(x)),dwMin) - - - newNfft = 2**nextpow2(np.ceil(wnNew/dwMin))+1 - if newNfft>nfft: - if (nfft<=2**15+1) and (newNfft>2**15+1): - warnings.warn('Spectrum matrix is very large (>33k). Memory problems may occur.') - - nfft = newNfft; - - - self.args = np.linspace(0,wnNew,nfft) - if method=='stineman': - self.data = stineman_interp(self.args,w,S1) - else: - intfun = interpolate.interp1d(w,S1,kind=method) - self.data = intfun(self.args) - self.data = self.data.clip(0) # clip negative values to 0 - - def normalize(self): - pass - def bandwidth(self,factors=0): - ''' Return some spectral bandwidth and irregularity factors - - - Returns - -------- - bw : arraylike - vector of bandwidth factors - - Parameters - ----------- - factors : array-like - Input vector 'factors' correspondence: - 0 alpha=m2/sqrt(m0*m4) (irregularity factor) - 1 eps2 = sqrt(m0*m2/m1^2-1) (narrowness factor) - 2 eps4 = sqrt(1-m2^2/(m0*m4))=sqrt(1-alpha^2) (broadness factor) - 3 Qp=(2/m0^2)int_0^inf f*S(f)^2 df (peakedness factor) - - Order of output is the same as order in 'factors' - - Example: - >>> import numpy as np - >>> import wafo.spectrum.models as sm - >>> Sj = sm.Jonswap(Hm0=3) - >>> w = np.linspace(0,4,256) - >>> S = SpecData1D(Sj(w),w) #Make spectrum object from numerical values - >>> S.bandwidth([0,1,2,3]) - array([ 0.65354446, 0.3975428 , 0.75688813, 2.00207912]) - ''' - - if self.freqtype in 'k': - vari = 'k'; - else: - vari = 'w'; - - m,mtxt = self.moment(nr=4,even=False) - sqrt = np.sqrt - fact = np.atleast_1d(factors) - alpha = m[2]/sqrt(m[0]*m[4]) - eps2 = sqrt(m[0]*m[2]/m[1]**2.-1.) - eps4 = sqrt(1.-m[2]**2./m[0]/m[4]); - f = self.args - S = self.data - Qp = 2/m[0]**2.*simps(f*S**2,x=f); - bw = np.array([alpha,eps2,eps4,Qp]) - return bw[fact] - - def characteristic(self,fact='Hm0',T=1200,g=9.81): - """Returns spectral characteristics and their covariance - - Parameters - ---------- - fact = vector with factor integers or a string or - a cellarray of strings, see below.(default [1]) - T : scalar - recording time (sec) (default 1200 sec = 20 min) - g : scalar - acceleration of gravity [m/s^2] - - Returns - ------- - ch = vector of spectral characteristics - R = matrix of the corresponding covariances given T - chtext = a list of strings describing the elements of ch, see example. - - - Description - ------------ - If input spectrum is of wave number type, output are factors for - corresponding 'k1D', else output are factors for 'freq'. - Input vector 'factors' correspondence: - 1 Hm0 = 4*sqrt(m0) Significant wave height - 2 Tm01 = 2*pi*m0/m1 Mean wave period - 3 Tm02 = 2*pi*sqrt(m0/m2) Mean zero-crossing period - 4 Tm24 = 2*pi*sqrt(m2/m4) Mean period between maxima - 5 Tm_10 = 2*pi*m_1/m0 Energy period - 6 Tp = 2*pi/{w | max(S(w))} Peak period - 7 Ss = 2*pi*Hm0/(g*Tm02^2) Significant wave steepness - 8 Sp = 2*pi*Hm0/(g*Tp^2) Average wave steepness - 9 Ka = abs(int S(w)*exp(i*w*Tm02) dw ) /m0 Groupiness parameter - 10 Rs = (S(0.092)+S(0.12)+S(0.15)/(3*max(S(w))) Quality control parameter - 11 Tp1 = 2*pi*int S(w)^4 dw Peak Period (robust estimate for Tp) - ------------------ - int w*S(w)^4 dw - - 12 alpha = m2/sqrt(m0*m4) Irregularity factor - 13 eps2 = sqrt(m0*m2/m1^2-1) Narrowness factor - 14 eps4 = sqrt(1-m2^2/(m0*m4))=sqrt(1-alpha^2) Broadness factor - 15 Qp = (2/m0^2)int_0^inf w*S(w)^2 dw Peakedness factor - - Order of output is same as order in 'factors' - The covariances are computed with a Taylor expansion technique - and is currently only available for factors 1, 2, and 3. Variances - are also available for factors 4,5,7,12,13,14 and 15 - - Quality control: - ---------------- - Critical value for quality control parameter Rs is Rscrit = 0.02 - for surface displacement records and Rscrit=0.0001 for records of - surface acceleration or slope. If Rs > Rscrit then probably there - are something wrong with the lower frequency part of S. - - Ss may be used as an indicator of major malfunction, by checking that - it is in the range of 1/20 to 1/16 which is the usual range for - locally generated wind seas. - - Examples: - --------- - >>> import numpy as np - >>> import wafo.spectrum.models as sm - >>> Sj = sm.Jonswap(Hm0=5) - >>> S = Sj.toSpecData() #Make spectrum ob - >>> S.characteristic(1) - (array([ 8.59007646]), array([[ 0.03040216]]), ['Tm01']) - - >>> [ch, R, txt] = S.characteristic([1,2,3]) # fact a vector of integers - >>> S.characteristic('Ss') # fact a string - (array([ 0.04963112]), array([[ 2.63624782e-06]]), ['Ss']) - - >>> S.characteristic(['Hm0','Tm02']) # fact a list of strings - (array([ 4.99833578, 8.03139757]), - array([[ 0.05292989, 0.02511371], - [ 0.02511371, 0.0274645 ]]), - ['Hm0', 'Tm02']) - - See also - --------- - bandwidth, - moment - - References - ---------- - Krogstad, H.E., Wolf, J., Thompson, S.P., and Wyatt, L.R. (1999) - 'Methods for intercomparison of wave measurements' - Coastal Enginering, Vol. 37, pp. 235--257 - - Krogstad, H.E. (1982) - 'On the covariance of the periodogram' - Journal of time series analysis, Vol. 3, No. 3, pp. 195--207 - - Tucker, M.J. (1993) - 'Recommended standard for wave data sampling and near-real-time processing' - Ocean Engineering, Vol.20, No.5, pp. 459--474 - - Young, I.R. (1999) - "Wind generated ocean waves" - Elsevier Ocean Engineering Book Series, Vol. 2, pp 239 - """ - - #% TODO % Need more checking on computing the variances for Tm24,alpha, eps2 and eps4 - #% TODO % Covariances between Tm24,alpha, eps2 and eps4 variables are also needed - NaN = np.nan - tfact = dict(Hm0=0,Tm01=1,Tm02=2,Tm24=3, Tm_10=4,Tp=5,Ss=6, Sp=7, Ka=8, - Rs=9, Tp1=10,Alpha=11,Eps2=12,Eps4=13,Qp=14) - tfact1 = ('Hm0','Tm01','Tm02','Tm24', 'Tm_10','Tp','Ss', 'Sp', 'Ka', - 'Rs', 'Tp1','Alpha','Eps2','Eps4','Qp') - - if isinstance(fact,str): - fact = list((fact,)) - if isinstance(fact,(list,tuple)): - nfact = [] - for k in fact: - if isinstance(k,str): - nfact.append(tfact.get(k.capitalize(),15)) - else: - nfact.append(k) - else: - nfact = fact; - - nfact = np.atleast_1d(nfact) - - if np.any((nfact>14) | (nfact<0)): - raise ValueError('Factor outside range (0,...,14)') - - vari = self.freqtype - - f = self.args.ravel() - S1 = self.data.ravel() - m,mtxt = self.moment(nr=4,even=False) - - #% moments corresponding to freq in Hz - for k in range(1,5): - m[k] = m[k]/(2*np.pi)**k - - pi = np.pi - ind = np.flatnonzero(f>0) - m.append(simps(S1[ind]/f[ind],f[ind])*2.*np.pi) #; % = m_1 - m_10 = simps(S1[ind]**2/f[ind],f[ind])*(2*pi)**2/T # % = COV(m_1,m0|T=t0) - m_11 = simps(S1[ind]**2./f[ind]**2,f[ind])*(2*pi)**3/T #% = COV(m_1,m_1|T=t0) - - sqrt = np.sqrt - #% Hm0 Tm01 Tm02 Tm24 Tm_10 - Hm0 = 4.*sqrt(m[0]); - Tm01 = m[0]/m[1]; - Tm02 = sqrt(m[0]/m[2]); - Tm24 = sqrt(m[2]/m[4]); - Tm_10= m[5]/m[0]; - - Tm12 = m[1]/m[2] - - ind = S1.argmax() - maxS = S1[ind] - #[maxS ind] = max(S1); - Tp = 2.*pi/f[ind] # % peak period /length - Ss = 2.*pi*Hm0/g/Tm02**2 # % Significant wave steepness - Sp = 2.*pi*Hm0/g/Tp**2 #; % Average wave steepness - Ka = abs(simps(S1*np.exp(1J*f*Tm02),f))/m[0]; #% groupiness factor - - #% Quality control parameter - #% critical value is approximately 0.02 for surface displacement records - #% If Rs>0.02 then there are something wrong with the lower frequency part - #% of S. - Rs = np.sum(np.interp(np.r_[0.0146, 0.0195, 0.0244]*2*pi,f,S1))/3./maxS; - Tp2 = 2*pi*simps(S1**4,f)/simps(f*S1**4,f); - - - alpha1 = Tm24/Tm02 #; % m(3)/sqrt(m(1)*m(5)); - eps2 = sqrt(Tm01/Tm12-1.)# % sqrt(m(1)*m(3)/m(2)^2-1); - eps4 = sqrt(1.-alpha1**2) #; % sqrt(1-m(3)^2/m(1)/m(5)); - Qp = 2./m[0]**2*simps(f*S1**2,f); - - - - - ch = np.r_[Hm0, Tm01, Tm02, Tm24, Tm_10, Tp, Ss, Sp, Ka, Rs, Tp2, alpha1, eps2, eps4, Qp] - - - #% Select the appropriate values - ch = ch[nfact] - chtxt = [tfact1[i] for i in nfact] - - #if nargout>1, - #% covariance between the moments: - #%COV(mi,mj |T=t0) = int f^(i+j)*S(f)^2 df/T - mij,mijtxt = self.moment(nr=8,even=False,j=1) - for ix,tmp in enumerate(mij): - mij[ix] = tmp/T/((2.*pi)**(ix-1.0)) - - - #% and the corresponding variances for - #%{'hm0', 'tm01', 'tm02', 'tm24', 'tm_10','tp','ss', 'sp', 'ka', 'rs', 'tp1','alpha','eps2','eps4','qp'} - R = np.r_[4*mij[0]/m[0], - mij[0]/m[1]**2.-2.*m[0]*mij[1]/m[1]**3.+m[0]**2.*mij[2]/m[1]**4., - 0.25*(mij[0]/(m[0]*m[2])-2.*mij[2]/m[2]**2+m[0]*mij[4]/m[2]**3), - 0.25*(mij[4]/(m[2]*m[4])-2*mij[6]/m[4]**2+m[2]*mij[8]/m[4]**3) , - m_11/m[0]**2+(m[5]/m[0]**2)**2*mij[0]-2*m[5]/m[0]**3*m_10, - NaN, - (8*pi/g)**2*(m[2]**2/(4*m[0]**3)*mij[0]+mij[4]/m[0]-m[2]/m[0]**2*mij[2]), - NaN*np.ones(4), - m[2]**2*mij[0]/(4*m[0]**3*m[4])+mij[4]/(m[0]*m[4])+mij[8]*m[2]**2/(4*m[0]*m[4]**3)- - m[2]*mij[2]/(m[0]**2*m[4])+m[2]**2*mij[4]/(2*m[0]**2*m[4]**2)-m[2]*mij[6]/m[0]/m[4]**2, - (m[2]**2*mij[0]/4+(m[0]*m[2]/m[1])**2*mij[2]+m[0]**2*mij[4]/4-m[2]**2*m[0]*mij[1]/m[1]+ - m[0]*m[2]*mij[2]/2-m[0]**2*m[2]/m[1]*mij[3])/eps2**2/m[1]**4, - (m[2]**2*mij[0]/(4*m[0]**2)+mij[4]+m[2]**2*mij[8]/(4*m[4]**2)-m[2]*mij[2]/m[0]+ - m[2]**2*mij[4]/(2*m[0]*m[4])-m[2]*mij[6]/m[4])*m[2]**2/(m[0]*m[4]*eps4)**2, - NaN]; - - #% and covariances by a taylor expansion technique: - #% Cov(Hm0,Tm01) Cov(Hm0,Tm02) Cov(Tm01,Tm02) - S0 = np.r_[ 2./(sqrt(m[0])*m[1])*(mij[0]-m[0]*mij[1]/m[1]), - 1./sqrt(m[2])*(mij[0]/m[0]-mij[2]/m[2]), - 1./(2*m[1])*sqrt(m[0]/m[2])*(mij[0]/m[0]-mij[2]/m[2]-mij[1]/m[1]+m[0]*mij[3]/(m[1]*m[2]))] - - R1 = np.ones((15,15)); - R1[:,:] = NaN - for ix,Ri in enumerate(R): - R1[ix,ix] = Ri - - - - R1[0,2:4] = S0[:2]; - R1[1,2] = S0[2]; - for ix in [0,1]: #%make lower triangular equal to upper triangular part - R1[ix+1:,ix] = R1[ix,ix+1:] - - - R = R[nfact] - R1= R1[nfact,:][:,nfact] - - - #% Needs further checking: - #% Var(Tm24)= 0.25*(mij[4]/(m[2]*m[4])-2*mij[6]/m[4]**2+m[2]*mij[8]/m[4]**3) ... - return ch, R1, chtxt - - def setlabels(self): - ''' Set automatic title, x-,y- and z- labels on SPECDATA object - - based on type, angletype, freqtype - ''' - - N = len(self.type); - if N==0: - raise ValueError('Object does not appear to be initialized, it is empty!') - - labels = ['','',''] - if self.type.endswith('dir'): - title = 'Directional Spectrum' - if self.freqtype.startswith('w'): - labels[0] = 'Frequency [rad/s]' - labels[2] = 'S(w,\theta) [m^2 s / rad^2]' - else: - labels[0] = 'Frequency [Hz]' - labels[2] = 'S(f,\theta) [m^2 s / rad]' - - if self.angletype.startswith('r'): - labels[1] = 'Wave directions [rad]' - elif self.angletype.startswith('d'): - labels[1] = 'Wave directions [deg]' - elif self.type.endswith('freq'): - title = 'Spectral density' - if self.freqtype.startswith('w'): - labels[0] = 'Frequency [rad/s]' - labels[1] = 'S(w) [m^2 s/ rad]' - else: - labels[0] = 'Frequency [Hz]' - labels[1] = 'S(f) [m^2 s]' - else: - title = 'Wave Number Spectrum' - labels[0] = 'Wave number [rad/m]' - if self.type.endswith('k1d'): - labels[1] = 'S(k) [m^3/ rad]' - elif self.type.endswith('k2d'): - labels[1] = labels[0] - labels[2] = 'S(k1,k2) [m^4/ rad^2]' - else: - raise ValueError('Object does not appear to be initialized, it is empty!') - if self.norm!=0: - title = 'Normalized ' + title - labels[0] = 'Normalized ' + labels[0].split('[')[0] - if not self.type.endswith('dir'): - labels[1] = labels[1].split('[')[0] - labels[2] = labels[2].split('[')[0] - - self.labels.title = title - self.labels.xlab = labels[0] - self.labels.ylab = labels[1] - self.labels.zlab = labels[2] -class SpecData2D(WafoData): - """ Container class for 2D spectrum data objects in WAFO - - Member variables - ---------------- - data : array_like - args : vector for 1D, list of vectors for 2D, 3D, ... - - type : string - spectrum type (default 'freq') - freqtype : letter - frequency type (default 'w') - angletype : string - angle type of directional spectrum (default 'radians') - - Examples - -------- - >>> import numpy as np - >>> import wafo.spectrum.models as sm - >>> Sj = sm.Jonswap(Hm0=3) - >>> w = np.linspace(0,4,256) - >>> S = SpecData1D(Sj(w),w) #Make spectrum object from numerical values - - See also - -------- - WafoData - CovData - """ - - def __init__(self,*args,**kwds): - super(SpecData2D, self).__init__(*args,**kwds) - - self.name='WAFO Spectrum Object' - self.type='freq' - self.freqtype='w' - self.angletype='' - self.h=np.inf - self.tr=None - self.phi=0. - self.v=0. - self.norm=0 - somekeys = ['angletype', 'phi', 'name', 'h', 'tr', 'freqtype', 'v', 'type', 'norm'] - - self.__dict__.update(sub_dict_select(kwds,somekeys)) - - if self.type.endswith('dir') and self.angletype=='': - self.angletype = 'radians' - - self.setlabels() - - def toacf(self): - pass - def sim(self): - pass - def sim_nl(self): - pass - def rotate(self): - pass - def moment(self,nr=2,vari='xt',even=True): - ''' Calculates spectral moments from spectrum - - Parameters - ---------- - nr : int - order of moments (maximum 4) - vari : string - variables in model, optional when two-dim.spectrum, - string with 'x' and/or 'y' and/or 't' - even : bool - False for all moments, - True for only even orders - - Returns - ------- - m : list of moments - mtext : list of strings describing the elements of m, see below - - Details - ------- - Calculates spectral moments of up to order four by use of - Simpson-integration. - - // - m_jkl=|| k1^j*k2^k*w^l S(w,th) dw dth - // - - where k1=w^2/gravity*cos(th-phi), k2=w^2/gravity*sin(th-phi) - and phi is the angle of the rotation in S.phi. If the spectrum - has field .g, gravity is replaced by S.g. - - The strings in output mtext have the same position in the cell array - as the corresponding numerical value has in output m - Notation in mtext: 'm0' is the variance, - 'mx' is the first-order moment in x, - 'mxx' is the second-order moment in x, - 'mxt' is the second-order cross moment between x and t, - 'myyyy' is the fourth-order moment in y - etc. - For the calculation of moments see Baxevani et al. - - Example: - S=demospec('dir') - [m,mtext]=spec2mom(S,2,'xyt') - - References - ---------- - Baxevani A. et al. (2001) - Velocities for Random Surfaces - ''' - -##% Tested on: Matlab 6.0 -##% Tested on: Matlab 5.3 -##% History: -##% Revised by I.R. 04.04.2001: Introducing the rotation angle phi. -##% Revised by A.B. 23.05.2001: Correcting 'mxxyy' and introducing -##% 'mxxyt','mxyyt' and 'mxytt'. -##% Revised by A.B. 21.10.2001: Correcting 'mxxyt'. -##% Revised by A.B. 21.10.2001: Adding odd-order moments. -##% By es 27.08.1999 - - - pi= np.pi - two_dim_spectra = ['dir','encdir','k2d'] - if self.type not in two_dim_spectra: - raise ValueError('Unknown 2D spectrum type!') - -## if (vari==None and nr<=1: -## vari='x' -## elif vari==None: -## vari='xt' -## else #% secure the mutual order ('xyt') -## vari=''.join(sorted(vari.lower())) -## Nv=len(vari) -## -## if vari[0]=='t' and Nv>1: -## vari = vari[1::]+ vari[0] -## -## Nv = len(vari) -## -## if not self.type.endswith('dir'): -## S1 = self.tospec(self.type[:-2]+'dir') -## else: -## S1 = self; -## w = np.ravel(S1.args[0]) -## theta = S1.args[1]-S1.phi -## S = S1.data -## Sw = simps(S,x=theta) -## m = [simps(Sw,x=w)] -## mtext=['m0'] -## -## if nr>0: -## -## nw=w.size -## if strcmpi(vari(1),'x') -## Sc=simpson(th,S1.S.*(cos(th)*ones(1,nw))).'; -## % integral S*cos(th) dth -## end -## if strcmpi(vari(1),'y') -## Ss=simpson(th,S1.S.*(sin(th)*ones(1,nw))).'; -## % integral S*sin(th) dth -## if strcmpi(vari(1),'x') -## Sc=simpson(th,S1.S.*(cos(th)*ones(1,nw))).'; -## end -## end -## if ~isfield(S1,'g') -## S1.g=gravity; -## end -## kx=w.^2/S1.g(1); % maybe different normalization in x and y => diff. g -## ky=w.^2/S1.g(end); -## -## if Nv>=1 -## switch vari -## case 'x' -## vec = kx.*Sc; -## mtext(end+1)={'mx'}; -## case 'y' -## vec = ky.*Ss; -## mtext(end+1)={'my'}; -## case 't' -## vec = w.*Sw; -## mtext(end+1)={'mt'}; -## end -## else -## vec = [kx.*Sc ky.*Ss w*Sw]; -## mtext(end+(1:3))={'mx', 'my', 'mt'}; -## end -## if nr>1 -## if strcmpi(vari(1),'x') -## Sc=simpson(th,S1.S.*(cos(th)*ones(1,nw))).'; -## % integral S*cos(th) dth -## Sc2=simpson(th,S1.S.*(cos(th).^2*ones(1,nw))).'; -## % integral S*cos(th)^2 dth -## end -## if strcmpi(vari(1),'y')||strcmpi(vari(2),'y') -## Ss=simpson(th,S1.S.*(sin(th)*ones(1,nw))).'; -## % integral S*sin(th) dth -## Ss2=simpson(th,S1.S.*(sin(th).^2*ones(1,nw))).'; -## % integral S*sin(th)^2 dth -## if strcmpi(vari(1),'x') -## Scs=simpson(th,S1.S.*((cos(th).*sin(th))*ones(1,nw))).'; -## % integral S*cos(th)*sin(th) dth -## end -## end -## if ~isfield(S1,'g') -## S1.g=gravity; -## end -## -## if Nv==2 -## switch vari -## case 'xy' -## vec=[kx.*Sc ky.*Ss kx.^2.*Sc2 ky.^2.*Ss2 kx.*ky.*Scs]; -## mtext(end+(1:5))={'mx','my','mxx', 'myy', 'mxy'}; -## case 'xt' -## vec=[kx.*Sc w.*Sw kx.^2.*Sc2 w.^2.*Sw kx.*w.*Sc]; -## mtext(end+(1:5))={'mx','mt','mxx', 'mtt', 'mxt'}; -## case 'yt' -## vec=[ky.*Ss w.*Sw ky.^2.*Ss2 w.^2.*Sw ky.*w.*Ss]; -## mtext(end+(1:5))={'my','mt','myy', 'mtt', 'myt'}; -## end -## else -## vec=[kx.*Sc ky.*Ss w.*Sw kx.^2.*Sc2 ky.^2.*Ss2 w.^2.*Sw kx.*ky.*Scs kx.*w.*Sc ky.*w.*Ss]; -## mtext(end+(1:9))={'mx','my','mt','mxx', 'myy', 'mtt', 'mxy', 'mxt', 'myt'}; -## end -## if nr>3 -## if strcmpi(vari(1),'x') -## Sc3=simpson(th,S1.S.*(cos(th).^3*ones(1,nw))).'; -## % integral S*cos(th)^3 dth -## Sc4=simpson(th,S1.S.*(cos(th).^4*ones(1,nw))).'; -## % integral S*cos(th)^4 dth -## end -## if strcmpi(vari(1),'y')||strcmpi(vari(2),'y') -## Ss3=simpson(th,S1.S.*(sin(th).^3*ones(1,nw))).'; -## % integral S*sin(th)^3 dth -## Ss4=simpson(th,S1.S.*(sin(th).^4*ones(1,nw))).'; -## % integral S*sin(th)^4 dth -## if strcmpi(vari(1),'x') %both x and y -## Sc2s=simpson(th,S1.S.*((cos(th).^2.*sin(th))*ones(1,nw))).'; -## % integral S*cos(th)^2*sin(th) dth -## Sc3s=simpson(th,S1.S.*((cos(th).^3.*sin(th))*ones(1,nw))).'; -## % integral S*cos(th)^3*sin(th) dth -## Scs2=simpson(th,S1.S.*((cos(th).*sin(th).^2)*ones(1,nw))).'; -## % integral S*cos(th)*sin(th)^2 dth -## Scs3=simpson(th,S1.S.*((cos(th).*sin(th).^3)*ones(1,nw))).'; -## % integral S*cos(th)*sin(th)^3 dth -## Sc2s2=simpson(th,S1.S.*((cos(th).^2.*sin(th).^2)*ones(1,nw))).'; -## % integral S*cos(th)^2*sin(th)^2 dth -## end -## end -## if Nv==2 -## switch vari -## case 'xy' -## vec=[vec kx.^4.*Sc4 ky.^4.*Ss4 kx.^3.*ky.*Sc3s ... -## kx.^2.*ky.^2.*Sc2s2 kx.*ky.^3.*Scs3]; -## mtext(end+(1:5))={'mxxxx','myyyy','mxxxy','mxxyy','mxyyy'}; -## case 'xt' -## vec=[vec kx.^4.*Sc4 w.^4.*Sw kx.^3.*w.*Sc3 ... -## kx.^2.*w.^2.*Sc2 kx.*w.^3.*Sc]; -## mtext(end+(1:5))={'mxxxx','mtttt','mxxxt','mxxtt','mxttt'}; -## case 'yt' -## vec=[vec ky.^4.*Ss4 w.^4.*Sw ky.^3.*w.*Ss3 ... -## ky.^2.*w.^2.*Ss2 ky.*w.^3.*Ss]; -## mtext(end+(1:5))={'myyyy','mtttt','myyyt','myytt','myttt'}; -## end -## else -## vec=[vec kx.^4.*Sc4 ky.^4.*Ss4 w.^4.*Sw kx.^3.*ky.*Sc3s ... -## kx.^2.*ky.^2.*Sc2s2 kx.*ky.^3.*Scs3 kx.^3.*w.*Sc3 ... -## kx.^2.*w.^2.*Sc2 kx.*w.^3.*Sc ky.^3.*w.*Ss3 ... -## ky.^2.*w.^2.*Ss2 ky.*w.^3.*Ss kx.^2.*ky.*w.*Sc2s ... -## kx.*ky.^2.*w.*Scs2 kx.*ky.*w.^2.*Scs]; -## mtext(end+(1:15))={'mxxxx','myyyy','mtttt','mxxxy','mxxyy',... -## 'mxyyy','mxxxt','mxxtt','mxttt','myyyt','myytt','myttt','mxxyt','mxyyt','mxytt'}; -## -## end % if Nv==2 ... else ... -## end % if nr>3 -## end % if nr>1 -## m=[m simpson(w,vec)]; -## end % if nr>0 -## % end; %%if Nv==1... else... to be removed -## end % ... else two-dim spectrum - - - - def interp(self): - pass - def normalize(self): - pass - def bandwidth(self): - pass - def setlabels(self): - ''' Set automatic title, x-,y- and z- labels on SPECDATA object - - based on type, angletype, freqtype - ''' - - N = len(self.type); - if N==0: - raise ValueError('Object does not appear to be initialized, it is empty!') - - labels = ['','',''] - if self.type.endswith('dir'): - title = 'Directional Spectrum' - if self.freqtype.startswith('w'): - labels[0] = 'Frequency [rad/s]' - labels[2] = 'S(w,\theta) [m**2 s / rad**2]' - else: - labels[0] = 'Frequency [Hz]' - labels[2] = 'S(f,\theta) [m**2 s / rad]' - - if self.angletype.startswith('r'): - labels[1] = 'Wave directions [rad]' - elif self.angletype.startswith('d'): - labels[1] = 'Wave directions [deg]' - elif self.type.endswith('freq'): - title = 'Spectral density' - if self.freqtype.startswith('w'): - labels[0] = 'Frequency [rad/s]' - labels[1] = 'S(w) [m**2 s/ rad]' - else: - labels[0] = 'Frequency [Hz]' - labels[1] = 'S(f) [m**2 s]' - else: - title = 'Wave Number Spectrum' - labels[0] = 'Wave number [rad/m]' - if self.type.endswith('k1d'): - labels[1] = 'S(k) [m**3/ rad]' - elif self.type.endswith('k2d'): - labels[1] = labels[0] - labels[2] = 'S(k1,k2) [m**4/ rad**2]' - else: - raise ValueError('Object does not appear to be initialized, it is empty!') - if self.norm!=0: - title = 'Normalized ' + title - labels[0] = 'Normalized ' + labels[0].split('[')[0] - if not self.type.endswith('dir'): - labels[1] = labels[1].split('[')[0] - labels[2] = labels[2].split('[')[0] - - self.labels.title = title - self.labels.xlab = labels[0] - self.labels.ylab = labels[1] - self.labels.zlab = labels[2] - -class CovData1D(WafoData): - """ Container class for 1D covariance data objects in WAFO - - Member variables - ---------------- - data : array_like - args : vector for 1D, list of vectors for 2D, 3D, ... - - type : string - spectrum type, one of 'freq', 'k1d', 'enc' (default 'freq') - lagtype : letter - lag type, one of: 'x', 'y' or 't' (default 't') - - - Examples - -------- - >>> import numpy as np - >>> import wafo.spectrum.models as sm - >>> Sj = sm.Jonswap(Hm0=3) - >>> w = np.linspace(0,4,256) - >>> S = SpecData1D(Sj(w),w) #Make spectrum object from numerical values - - See also - -------- - WafoData - CovData - """ - - def __init__(self,*args,**kwds): - super(CovData1D, self).__init__(*args,**kwds) - - self.name='WAFO Covariance Object' - self.type='time' - self.lagtype='t' - self.h=np.inf - self.tr=None - self.phi=0. - self.v=0. - self.norm=0 - somekeys = ['phi', 'name', 'h', 'tr', 'lagtype', 'v', 'type', 'norm'] - - self.__dict__.update(sub_dict_select(kwds,somekeys)) - - #self.setlabels() - def copy(self): - kwds = self.__dict__.copy() - wdata = CovData1D(**kwds) - return wdata - - def tospec(self,rate=None,method='linear',nugget=0.0,trunc=1e-5,fast=True): - '''Computes spectral density from the auto covariance function - - Parameters - ---------- - rate = scalar, int - 1,2,4,8...2^r, interpolation rate for f (default 1) - - method: string - interpolation method 'stineman', 'linear', 'cubic' - - nugget = scalar, real - nugget effect to ensure that round off errors do not result in - negative spectral estimates. Good choice might be 10^-12. - - trunc : scalar, real - truncates all spectral values where S/max(S) < trunc - 0 <= trunc <1 This is to ensure that high frequency - noise is not added to the spectrum. (default 1e-5) - fast : bool - if True : zero-pad to obtain power of 2 length ACF (default) - otherwise no zero-padding of ACF, slower but more accurate. - - Returns - -------- - S = SpecData1D object - spectral density - - NB! This routine requires that the covariance is evenly spaced - starting from zero lag. Currently only capable of 1D matrices. - - Example: - >>> import wafo.spectrum.models as sm - >>> import numpy as np - >>> import scipy.signal.signaltools as st - >>> L = 129 - >>> t = np.linspace(0,75,L) - >>> R = np.zeros(L) - >>> win = st.parzen(41) - >>> R[0:20] = win[20:41] - >>> R0 = CovData1D(R,t) - >>> S0 = R0.tospec() - - >>> Sj = sm.Jonswap() - >>> S = Sj.toSpecData() - >>> R2 = S.toacf() - >>> S1 = R2.tospec() - >>> assert(all(abs(S1.data-S.data)<1e-4) ,'COV2SPEC') - - See also - -------- - spec2cov - datastructures - ''' - - dT = self.sampling_period() - # dT = time-step between data points.(default = T(2)-T1)). - - ACF, ti = np.atleast_1d(self.data,self.args) - - if self.lagtype in 't': - spectype = 'freq' - ftype = 'w' - else: - spectype = 'k1d' - ftype = 'k' - - if rate is None: - rate = 1 #;%interpolation rate - else: - rate = 2**nextpow2(rate) #;%make sure rate is a power of 2 - - - #% add a nugget effect to ensure that round off errors - #% do not result in negative spectral estimates - ACF[0] = ACF[0] +nugget - n = ACF.size - # embedding a circulant vector and Fourier transform - if fast: - nfft = 2**nextpow2(2*n-2) - else: - nfft = 2*n-2 - - nf = nfft/2 #;% number of frequencies - fft = np.fft.fft - ACF = np.r_[ACF,np.zeros(nfft-2*n+2),ACF[n-1:0:-1]] - - Rper = (fft(ACF,nfft).real).clip(0) #% periodogram - RperMax = Rper.max() - Rper = np.where(Rper1: - So.args = np.linspace(0,pi/dT,nf*rate); - if method=='stineman': - So.data = stineman_interp(So.args,w,S) - else: - intfun = interpolate.interp1d(w,S,kind=method) - So.data = intfun(So.args) - So.data = So.data.clip(0) # clip negative values to 0 - - return So - - - - - - def sampling_period(self): - ''' Returns sampling interval - - Returns - --------- - dT : scalar - sampling interval, unit: - [s] if lagtype=='t' - [m] otherwise - - - See also - ''' - dt1 = self.args[1]-self.args[0] - N = np.size(self.args)-1 - T = self.args[-1]-self.args[0] - dt = T/N - - return dt - - def sim(self,ns=None,cases=1,dt=None,iseed=None,derivative=False): - ''' Simulates a Gaussian process and its derivative from ACF - - Parameters - ---------- - ns : scalar - number of simulated points. (default length(S)-1=n-1). - If ns>n-1 it is assummed that R(k)=0 for all k>n-1 - cases : scalar - number of replicates (default=1) - dt : scalar - step in grid (default dt is defined by the Nyquist freq) - iseed : int or state - starting state/seed number for the random number generator - (default none is set) - derivative : bool - if true : return derivative of simulated signal as well - otherwise - - - Returns - ------- - xs = a cases+1 column matrix ( t,X1(t) X2(t) ...). - xsder = a cases+1 column matrix ( t,X1'(t) X2'(t) ...). - - Details - ------- - Performs a fast and exact simulation of stationary zero mean - Gaussian process through circulant embedding of the covariance matrix. - - If the ACF has a non-empty field .tr, then the transformation is - applied to the simulated data, the result is a simulation of a transformed - Gaussian process. - - Note: The simulation may give high frequency ripple when used with a - small dt. - - Example: - >>> import wafo.spectrum.models as sm - >>> Sj = sm.Jonswap() - >>> S = Sj.toSpecData() #Make spec - >>> R = S.toacf() - >>> x = R.sim(ns=1000,dt=0.2) - - See also - -------- - spec2sdat, gaus2dat - - Reference - ----------- - C.R Dietrich and G. N. Newsam (1997) - "Fast and exact simulation of stationary - Gaussian process through circulant embedding - of the Covariance matrix" - SIAM J. SCI. COMPT. Vol 18, No 4, pp. 1088-1107 - - ''' - - # TODO fix it, it does not work - #% add a nugget effect to ensure that round off errors - #% do not result in negative spectral estimates - nugget=0 #;%10^-12; - - _set_seed(iseed) - - ACF = self.data.ravel() - n = ACF.size - - - I = ACF.argmax() - if I!=0: - raise ValueError('ACF does not have a maximum at zero lag') - - ACF.shape = (n,1) - - dT = self.sampling_period() - - - fft = np.fft.fft - - x=np.zeros((ns,cases+1)) - - - if derivative: - xder=x.copy() - - - #% add a nugget effect to ensure that round off errors - #% do not result in negative spectral estimates - ACF[0]=ACF[0]+nugget; - - #% Fast and exact simulation of simulation of stationary - #% Gaussian process throug circulant embedding of the - #% Covariance matrix - floatinfo = numpy.finfo(float) - if (abs(ACF[-1])>floatinfo.eps): #% assuming ACF(n+1)==0 - m2=2*n-1; - nfft=2**nextpow2(max(m2,2*ns)); - ACF=np.r_[ACF,np.zeros((nfft-m2,1)),ACF[-1:0:-1,:]] - - #disp('Warning: I am now assuming that ACF(k)=0 ') - #disp('for k>MAXLAG.') - - else: # % ACF(n)==0 - m2=2*n-2; - nfft=2**nextpow2(max(m2,2*ns)) - ACF=np.r_[ACF,np.zeros((nfft-m2,1)),ACF[n-1:1:-1,:]] - - #%m2=2*n-2; - - S=fft(ACF,nfft,axis=0).real #;% periodogram - - - I=S.argmax() - k=np.flatnonzero(S<0); - if k.size>0: - #disp('Warning: Not able to construct a nonnegative circulant ') - #disp('vector from the ACF. Apply the parzen windowfunction ') - #disp('to the ACF in order to avoid this.') - #disp('The returned result is now only an approximation.') - - # truncating negative values to zero to ensure that - # that this noise is not added to the simulated timeseries - - S[k] = 0. - - ix = np.flatnonzero(k>2*I) - if ix.size>0: -## % truncating all oscillating values above 2 times the peak -## % frequency to zero to ensure that -## % that high frequency noise is not added to -## % the simulated timeseries. - ix0 = k[ix[0]] - S[ix0:-ix0] =0.0 - - - sqrt = np.sqrt - trunc = 1e-5; - maxS = S[I] - k=np.flatnonzero(S[I:-I]0: - S[k+I]=0. - #% truncating small values to zero to ensure that - #% that high frequency noise is not added to - #% the simulated timeseries - - cases1 = np.floor(cases/2); - cases2 = np.ceil(cases/2); -#% Generate standard normal random numbers for the simulations -#% ----------------------------------------------------------- - randn = np.random.randn - epsi = randn(nfft,cases2)+1j*randn(nfft,cases2); - Ssqr=sqrt(S/(nfft)) #; %sqrt(S(wn)*dw ) - ephat=epsi*Ssqr[:,np.newaxis] - y=fft(ephat,nfft); - x[:,1:cases+1]=np.hstack((y[2:ns+2,0:cases2].real, y[2:ns+2,0:cases1].imag)) - - - - x[:,0]=np.linspace(0,(ns-1)*dT,ns) #%(0:dT:(dT*(np-1)))'; - - if derivative: - Ssqr = Ssqr*np.r_[0:(nfft/2+1),-(nfft/2-1):0]*2*np.pi/nfft/dT - ephat = epsi*Ssqr[:,np.newaxis] - y = fft(ephat,nfft); - xder[:,1:(cases+1)]= np.hstack((y[2:ns+2,0:cases2].imag -y[2:ns+2,0:cases1].real)) - xder[:,0]=x[:,0] - - if self.tr is not None: - np.disp(' Transforming data.') - g=self.tr; - G=np.fliplr(g); #% the invers of g - if derivative: - for ix in range(cases): - tmp=tranproc(np.hstack((x[:,ix+1], xder[:,ix+1])),G); - x[:,ix+1]=tmp[:,0]; - xder[:,ix+1]=tmp[:,1]; - - else: - for ix in range(cases): - x[:,ix+1]=tranproc(x[:,ix+1],G); - - if derivative: - return x,xder - else: - return x - - - -class TimeSeries(WafoData): - ''' Container class for 1D TimeSeries data objects in WAFO - - Member variables - ---------------- - data : array_like - args : vector for 1D, list of vectors for 2D, 3D, ... - - type : integer or string - sensor type for time series (default 'n' : Surface elevation) - see sensortype for more options - position : vector of size 3 - instrument position relative to the coordinate system - - ''' - def __init__(self,*args,**kwds): - super(TimeSeries, self).__init__(*args,**kwds) - self.name='WAFO TimeSeries Object' - self.type='n' - self.position = np.zeros(3) - somekeys = ['type', 'position'] - self.__dict__.update(sub_dict_select(kwds,somekeys)) - - #self.setlabels() - if not any(self.args): - n = len(self.data) - self.args = xrange(0,n) - - def sampling_period(self): - ''' Returns sampling interval - - Returns - --------- - dT : scalar - sampling interval, unit: - [s] if lagtype=='t' - [m] otherwise - - - See also - ''' - dt1 = self.args[1]-self.args[0] - N = np.size(self.args)-1 - T = self.args[-1]-self.args[0] - dt = T/N - - return dt - - def acf(self,lag=None,flag='biased',norm=False,dt = None): - ''' Return auto covariance function from data. - - - R = ACF vector length L+1 - t = time lags length L+1 - stdev = estimated large lag standard deviation of the estimate - assuming x is a Gaussian process: - if R(k)=0 for all lags k>q then an approximation - of the variance for large samples due to Bartlett - var(R(k))=1/N*(R(0)^2+2*R(1)^2+2*R(2)^2+ ..+2*R(q)^2) - for k>q and where N=length(x). Special case is - white noise where it equals R(0)^2/N for k>0 - norm = 0 indicating that R is not normalized - - x = a column data vector or two column data matrix with sampled times and values. - L = the maximum time-lag for which the ACF is estimated. - (Default L=n-1) - plotflag = 1 then the ACF is plotted vs lag - 2 then the ACF is plotted vs lag in seconds - 3 then the ACF is plotted vs lag and vs lag (sec - dT = time-step between data points (default xn(2,1)-xn(1,1) or 1 Hz). - flag = 'biased' : scales the raw cross-correlation by 1/n. (default) - 'unbiased': scales the raw correlation by 1/(n-abs(k)), - where k is the index into the result. - - - - x may contain NaN's (i.e. missing values). - - Example: - x = load('sea.dat'); - rf = dat2cov(x,150,2) - ''' - n = len(self.data) - if not lag: - lag = n-1 - - x = self.data.copy() - indnan = numpy.isnan(x) - if any(indnan): - x = x - numpy.mean(x[1-indnan]) # remove the mean pab 09.10.2000 - #indnan = find(indnan); - Ncens = n - sum(indnan) - x[indnan] = 0. # pab 09.10.2000 much faster for censored samples - else: - indnan = None - Ncens = n - x = x - numpy.mean(x) - - fft = numpy.fft.fft - nfft = 2**nextpow2(n) - Rper = abs(fft(x,nfft))**2/Ncens # Raw periodogram - - R = numpy.real(fft(Rper))/nfft # %ifft=fft/nfft since Rper is real! - lags = np.arange(0,lag+1) - if flag.startswith('unbiased'): - # unbiased result, i.e. divide by n-abs(lag) - R=R[lags]*Ncens/ np.arange(Ncens,Ncens-lag,-1) - #else % biased result, i.e. divide by n - # r=r(1:L+1)*Ncens/Ncens; - - #c0 = R[0] - dT = self.sampling_period() - t = numpy.linspace(0,lag*dT,lag+1) - cumsum = numpy.cumsum - acf = CovData1D(R[lags],t) - acf.stdev=numpy.sqrt(numpy.r_[ 0, 1 ,1+2*cumsum(R[1:]**2)]/Ncens) - acf.children = [WafoData(-2.*acf.stdev[lags],t),WafoData(2.*acf.stdev[lags],t)] - return acf - - def spec(self): - pass - def turning_points(self,h=0,wavetype=None,output='tp'): - ''' Return turning points (tp) from data, optionally rainflowfiltered. - - Parameters - ---------- - h : scalar - a threshold; - if h<0, then tp=x; - if h=0, then tp is a sequence of turning points (default); - if h>0, then all rainflow cycles with height smaller than - h are removed. - - wavetype : string - defines the type of wave. Possible options are - 'mw' 'Mw' or 'none'. - If None all rainflow filtered min and max - will be returned, otherwise only the rainflow filtered - min and max, which define a wave according to the - wave definition, will be returned. - output : string - 'tp' if only returning tp - 'all' if tp and ind - - Returns - ------- - tp : TimeSeries object - with times and turning points. - - ind : arraylike - indices to the turning points in the original sequence. - - - Example: - x = load('sea.dat'); x1 = x(1:200,:); - tp = dat2tp(x1,0,'Mw'); tph = dat2tp(x1,0.3,'Mw'); - plot(x1(:,1),x1(:,2),tp(:,1),tp(:,2),'ro',tph(:,1),tph(:,2),'k*') - - See also - --------- - findcross, findrfc - ''' - - if h<0: - tp = self.copy() - ind = np.arange(tp.args.size) - if output.startswith('tp'): - return tp - else: - return tp, ind - x2 = self.data - - n = len(x2) - ind = findextrema(x2) - - if ind.size<2: - tp = None - ind = None - return tp - - - #% In order to get the exact up-crossing intensity from rfc by - #% mm2lc(tp2mm(rfc)) we have to add the indices - #% to the last value (and also the first if the - #% sequence of turning points does not start with a minimum). - - if x2[ind[0]]>x2[ind[1]]: - #% adds indices to first and last value - ind = np.r_[0, ind ,n] - else: # adds index to the last value - ind = np.r_[ind, n] - - - if h>0: - ind1 = findrfc(x2[ind],h) - ind = ind[ind1] - - Nm =ind.size #% number of min and Max - - if wavetype in ('mw','Mw'): - - xor = lambda a,b : a^b - #% make sure that the first is a Max if wdef == 'Mw' - #% or make sure that the first is a min if wdef == 'mw' - removeFirst = xor((x2[ind[0]]>x2[ind[1]]),wavetype.startswith('Mw')) - if removeFirst: - ind = ind[1::] - Nm = Nm-1; - - - #% make sure the number of minima and Maxima are according to the wavedef. - #% i.e., make sure Nm=length(ind) is odd - if (mod(Nm,2))!=1: - ind = ind[:-1] - Nm = Nm-1 - try: - t = self.args - except: - t = ind - tp = TimeSeries(self.data[ind],t) - - return tp - - def lc_spectrum(self): - pass - def cycles(self): - pass - def trough_crest(self): - pass - def wave_period(self): - pass - def reconstruct(self): - pass - def plot_wave_idx(self): - ''' spwaveplot - ''' - pass - def findoutliers(self): - pass - - -def sensortypeid(*sensortypes): - ''' Return ID for sensortype name - - Parameter - --------- - sensortypes : list of strings defining the sensortype - - Returns - ------- - sensorids : list of integers defining the sensortype - - Valid senor-ids and -types for time series are as follows: - 0, 'n' : Surface elevation (n=Eta) - 1, 'n_t' : Vertical surface velocity - 2, 'n_tt' : Vertical surface acceleration - 3, 'n_x' : Surface slope in x-direction - 4, 'n_y' : Surface slope in y-direction - 5, 'n_xx' : Surface curvature in x-direction - 6, 'n_yy' : Surface curvature in y-direction - 7, 'n_xy' : Surface curvature in xy-direction - 8, 'P' : Pressure fluctuation about static MWL pressure - 9, 'U' : Water particle velocity in x-direction - 10, 'V' : Water particle velocity in y-direction - 11, 'W' : Water particle velocity in z-direction - 12, 'U_t' : Water particle acceleration in x-direction - 13, 'V_t' : Water particle acceleration in y-direction - 14, 'W_t' : Water particle acceleration in z-direction - 15, 'X_p' : Water particle displacement in x-direction from its mean position - 16, 'Y_p' : Water particle displacement in y-direction from its mean position - 17, 'Z_p' : Water particle displacement in z-direction from its mean position - - Example: - >>> sensortypeid('W','v') - [11, 10] - >>> sensortypeid('rubbish') - [-1.#IND] - - See also sensortype, id - ''' - - sensorid_table = dict(n=0,n_t=1,n_tt=2,n_x=3,n_y=4,n_xx=5, - n_yy=6,n_xy=7,p=8,u=9,v=10,w=11,u_t=12, - v_t=13,w_t=14,x_p=15,y_p=16,z_p=17) - try: - return [sensorid_table.get(name.lower(),np.NAN) for name in sensortypes] - except: - raise ValueError('Input must be a string!') - - - -def sensortype(*sensorids): - ''' Return sensortype name - - Parameter - --------- - sensorids : vector or list of integers defining the sensortype - - Returns - ------- - sensornames : list of strings defining the sensortype - - Valid senor-ids and -types for time series are as follows: - 0, 'n' : Surface elevation (n=Eta) - 1, 'n_t' : Vertical surface velocity - 2, 'n_tt' : Vertical surface acceleration - 3, 'n_x' : Surface slope in x-direction - 4, 'n_y' : Surface slope in y-direction - 5, 'n_xx' : Surface curvature in x-direction - 6, 'n_yy' : Surface curvature in y-direction - 7, 'n_xy' : Surface curvature in xy-direction - 8, 'P' : Pressure fluctuation about static MWL pressure - 9, 'U' : Water particle velocity in x-direction - 10, 'V' : Water particle velocity in y-direction - 11, 'W' : Water particle velocity in z-direction - 12, 'U_t' : Water particle acceleration in x-direction - 13, 'V_t' : Water particle acceleration in y-direction - 14, 'W_t' : Water particle acceleration in z-direction - 15, 'X_p' : Water particle displacement in x-direction from its mean position - 16, 'Y_p' : Water particle displacement in y-direction from its mean position - 17, 'Z_p' : Water particle displacement in z-direction from its mean position - - Example: - >>> sensortype(range(3)) - ['n', 'n_t', 'n_tt'] - - See also sensortypeid, tran - ''' - - - - validNames = ('n','n_t','n_tt','n_x','n_y','n_xx', - 'n_yy','n_xy','p','u','v','w','u_t', - 'v_t','w_t','x_p','y_p','z_p',np.NAN) - ids = np.atleast_1d(*sensorids) - if isinstance(ids,list): - ids = np.hstack(ids) - N = len(validNames)-1 - - try: - return [validNames[id] for id in ids.clip(0,N)] - except: - raise ValueError('Input must be an integer!') - - -def main(): - from wafo.spectrum import models as sm - sensortype(range(21)) - w = np.linspace(0,3,100) - Sj = sm.Jonswap() - S = Sj.toSpecData() - #S = SpecData1D(Sj(w),w) - R = S.toacf(nr=1) - S1 = S.copy() - Si = R.tospec() - ns =5000; - dt = .2; - x1 = S.sim_nl(ns= ns,dt=dt); - x2 = TimeSeries(x1[:,1],x1[:,0]) - R = x2.acf(lag=100) - R.plot() - - S.plot('ro') - t = S.moment() - t1 = S.bandwidth([0,1,2,3]) - S1 = S.copy() - S1.resample(dt=0.3,method='cubic') - S1.plot('k+') - x = S1.sim(ns=100) - import pylab - pylab.clf() - pylab.plot(x[:,0],x[:,1]) - pylab.show() - - pylab.close('all') - np.disp('done') - - -if __name__ == '__main__': - if False: #True: # - import doctest - doctest.testmod() - else: - #main() - import wafo.spectrum.models as sm - Sj = sm.Jonswap(); - S = Sj.toSpecData(); - - R = S.toacf() - x = R.sim(ns=1000,dt=0.2) - S.characteristic(['hm0','tm02']) - ns =1000; dt = .2; - x1 = S.sim_nl(ns,dt=dt); - - x = np.arange(-2,2,0.2) - - # Plot 2 objects in one call - d2 = WafoData(np.sin(x),x,xlab='x',ylab='sin',title='sinus') - - - d0 = d2.copy() - d0.data = d0.data*0.9 - d1 = d2.copy() - d1.data = d1.data*1.2 - d1.children = [d0] - d2.children = [d1] - - d2.plot() - print 'Done' \ No newline at end of file diff --git a/pywafo/src/wafo/integrate.py b/pywafo/src/wafo/integrate.py index 3d88b96..4ca3ac3 100644 --- a/pywafo/src/wafo/integrate.py +++ b/pywafo/src/wafo/integrate.py @@ -6,9 +6,8 @@ from numpy import pi, sqrt, ones, zeros #@UnresolvedImport from scipy import integrate as intg import scipy.special.orthogonal as ort from scipy import special as sp -import matplotlib import pylab as plb -matplotlib.interactive(True) + _POINTS_AND_WEIGHTS = {} def humps(x=None): @@ -344,7 +343,7 @@ def h_roots(n, method='newton'): >>> import numpy as np >>> [x,w] = h_roots(10) >>> np.sum(x*w) - -5.2516042729766621e-019 + -5.2516042729766621e-19 See also -------- @@ -453,7 +452,7 @@ def j_roots(n, alpha, beta, method='newton'): -------- >>> [x,w]= j_roots(10,0,0) >>> sum(x*w) - 2.7755575615628914e-016 + 2.7755575615628914e-16 See also -------- @@ -554,7 +553,7 @@ def la_roots(n, alpha=0, method='newton'): >>> import numpy as np >>> [x,w] = h_roots(10) >>> np.sum(x*w) - -5.2516042729766621e-019 + -5.2516042729766621e-19 See also -------- @@ -939,7 +938,7 @@ def gaussq(fun, a, b, reltol=1e-3, abstol=1e-3, alpha=0, beta=0, wfun=1, integration of x**2 from 0 to 2 and from 1 to 4 >>> from scitools import numpyutils as npt - scitools.easyviz backend is gnuplot + scitools.easyviz backend is matplotlib >>> A = [0, 1]; B = [2,4] >>> fun = npt.wrap2callable('x**2') >>> [val1,err1] = gaussq(fun,A,B) @@ -1165,22 +1164,22 @@ def quadgr(fun,a,b,abseps=1e-5): >>> import numpy as np >>> Q, err = quadgr(np.log,0,1) >>> quadgr(np.exp,0,9999*1j*np.pi) - (-2.0000000000122617, 2.1933275196062141e-009) + (-2.0000000000122617, 2.1933275196062141e-09) >>> quadgr(lambda x: np.sqrt(4-x**2),0,2,1e-12) - (3.1415926535897811, 1.5809575870662229e-013) + (3.1415926535897811, 1.5809575870662229e-13) >>> quadgr(lambda x: x**-0.75,0,1) - (4.0000000000000266, 5.6843418860808015e-014) + (4.0000000000000266, 5.6843418860808015e-14) >>> quadgr(lambda x: 1./np.sqrt(1-x**2),-1,1) - (3.141596056985029, 6.2146261559092864e-006) + (3.141596056985029, 6.2146261559092864e-06) >>> quadgr(lambda x: np.exp(-x**2),-np.inf,np.inf,1e-9) #% sqrt(pi) - (1.7724538509055152, 1.9722334876348668e-011) + (1.7724538509055152, 1.9722334876348668e-11) >>> quadgr(lambda x: np.cos(x)*np.exp(-x),0,np.inf,1e-9) - (0.50000000000000044, 7.3296813063450372e-011) + (0.50000000000000044, 7.3296813063450372e-11) See also -------- @@ -1418,7 +1417,7 @@ def qdemo(f,a,b): #[x, w]=qrule(n,1) #x = (b-a)/2*x + (a+b)/2 % Transform base points X. #w = (b-a)/2*w % Adjust weigths. - #q = sum(feval(f,x).*w) + #q = sum(feval(f,x)*w) qg[k] = q eg[k] = abs(q - true_val) diff --git a/pywafo/src/wafo/kdetools.py b/pywafo/src/wafo/kdetools.py index 9a7f96c..fb9acd9 100644 --- a/pywafo/src/wafo/kdetools.py +++ b/pywafo/src/wafo/kdetools.py @@ -15,7 +15,7 @@ from numpy import pi, atleast_2d #@UnresolvedImport from misc import tranproc, trangood def sphere_volume(d, r=1.0): """ - Returns volume of d-dimensional sphere with radius r + Returns volume of d-dimensional sphere with radius r Parameters ---------- diff --git a/pywafo/src/wafo/objects.py b/pywafo/src/wafo/objects.py index 77f1ad1..dda6761 100644 --- a/pywafo/src/wafo/objects.py +++ b/pywafo/src/wafo/objects.py @@ -694,8 +694,9 @@ class TimeSeries(WafoData): Examples -------- >>> import wafo.data + >>> import wafo.objects as wo >>> x = wafo.data.sea() - >>> ts = wafo.objects.mat2timeseries(x) + >>> ts = wo.mat2timeseries(x) >>> rf = ts.tocovdata(lag=150) >>> h = rf.plot() @@ -771,8 +772,9 @@ class TimeSeries(WafoData): Example: -------- >>> import wafo.data + >>> import wafo.objects as wo >>> x = wafo.data.sea() - >>> ts = mat2timeseries(x) + >>> ts = wo.mat2timeseries(x) >>> acf = ts.tocovdata(150) >>> h = acf.plot() ''' @@ -814,6 +816,7 @@ class TimeSeries(WafoData): acf = _wafocov.CovData1D(R[lags], t) acf.stdev = sqrt(r_[ 0, 1 , 1 + 2 * cumsum(R[1:] ** 2)] / Ncens) acf.children = [WafoData(-2. * acf.stdev[lags], t), WafoData(2. * acf.stdev[lags], t)] + acf.plot_args_children = ['r:'] acf.norm = norm return acf diff --git a/pywafo/src/wafo/spectrum/__init__.py b/pywafo/src/wafo/spectrum/__init__.py index 9cfd44c..e56bd4c 100644 --- a/pywafo/src/wafo/spectrum/__init__.py +++ b/pywafo/src/wafo/spectrum/__init__.py @@ -8,4 +8,3 @@ Spectrum package in WAFO Toolbox. from wafo.spectrum.core import SpecData1D import wafo.spectrum.models import wafo.spectrum.dispersion_relation -#from wafo.data_structures import SpecData1D \ No newline at end of file diff --git a/pywafo/src/wafo/spectrum/core.py b/pywafo/src/wafo/spectrum/core.py index 30a6bb7..1678fdb 100644 --- a/pywafo/src/wafo/spectrum/core.py +++ b/pywafo/src/wafo/spectrum/core.py @@ -1,9 +1,9 @@ from __future__ import division -from scipy.misc.ppimport import ppimport +from wafo.misc import meshgrid from wafo.objects import mat2timeseries, TimeSeries import warnings import numpy as np -from numpy import (pi, inf, meshgrid, zeros, ones, where, nonzero, #@UnresolvedImport +from numpy import (pi, inf, zeros, ones, where, nonzero, #@UnresolvedImport flatnonzero, ceil, sqrt, exp, log, arctan2, #@UnresolvedImport tanh, cosh, sinh, random, atleast_1d, maximum, #@UnresolvedImport minimum, diff, isnan, any, r_, conj, mod, #@UnresolvedImport @@ -945,7 +945,8 @@ class SpecData1D(WafoData): ... res = fun(x2[:,1::],axis=0) ... m = res.mean() ... sa = res.std() - ... assert(np.abs(m-trueval)3, ## %compute the sum and frequency effects separately ## [svec, dvec] = disufq((amp.'),w,kw,min(h,10^30),g,nmin,nmax) @@ -1402,7 +1401,10 @@ class SpecData1D(WafoData): See also --------- - hermitetr, ochitr, lc2tr, dat2tr + transform.TrHermite + transform.TrOchi + objects.LevelCrossings.trdata + objects.TimeSeries.trdata References: ----------- @@ -1429,7 +1431,7 @@ class SpecData1D(WafoData): w = ravel(self.args) S = ravel(self.data) if self.freqtype in ['f', 'w']: - vari = 't' + #vari = 't' if self.freqtype == 'f': w = 2. * pi * w S = S / (2. * pi) @@ -1514,7 +1516,7 @@ class SpecData1D(WafoData): ## skew = sum((6*C2+8*E2).*E)/sa^3 % skewness ## kurt = 3+48*sum((C2+E2).*E2)/sa^4 % kurtosis return output - def testgaussian(self, ns,test0=None, cases=100, method='nonlinear',**opt): + def testgaussian(self, ns,test0=None, cases=100, method='nonlinear',verbose=False,**opt): ''' TESTGAUSSIAN Test if a stochastic process is Gaussian. @@ -1607,7 +1609,8 @@ class SpecData1D(WafoData): #xs = cov2sdat(R,[ns Nstep]); #[g, tmp] = dat2tr(xs,method, **opt); #test1 = [test1; tmp(:)] - print('finished %d of %d ' % (ix+1,rep) ) + if verbose: + print('finished %d of %d ' % (ix+1,rep) ) if rep>1: xs = acf.sim(ns=ns, cases=np.remainder(cases,rep)) @@ -1704,6 +1707,14 @@ class SpecData1D(WafoData): def nyquist_freq(self): """ Return Nyquist frequency + + Example + ------- + >>> import wafo.spectrum.models as sm + >>> Sj = sm.Jonswap(Hm0=5) + >>> S = Sj.tospecdata() #Make spectrum ob + >>> S.nyquist_freq() + 3.0 """ return self.args[-1] @@ -1722,8 +1733,11 @@ class SpecData1D(WafoData): Example ------- - S = jonswap - dt = spec2dt(S) + >>> import wafo.spectrum.models as sm + >>> Sj = sm.Jonswap(Hm0=5) + >>> S = Sj.tospecdata() #Make spectrum ob + >>> S.sampling_period() + 1.0471975511965976 See also ''' @@ -1880,9 +1894,16 @@ class SpecData1D(WafoData): Example: ------- - S = jonswap - [Sn,mn4] = specnorm(S) - mts = spec2mom(S,2) % Should be equal to one! + >>> import wafo.spectrum.models as sm + >>> Sj = sm.Jonswap(Hm0=5) + >>> S = Sj.tospecdata() #Make spectrum ob + >>> S.moment(2) + ([1.5614600345079888, 0.95567089481941048], ['m0', 'm0tt']) + >>> Sn = S.copy(); Sn.normalize() + + Now the moments should be one + >>> Sn.moment(2) + ([1.0000000000000004, 0.99999999999999967], ['m0', 'm0tt']) ''' mom, unused_mtext = self.moment(nr=4, even=True) m0 = mom[0] @@ -2018,7 +2039,6 @@ class SpecData1D(WafoData): Examples: --------- - >>> import numpy as np >>> import wafo.spectrum.models as sm >>> Sj = sm.Jonswap(Hm0=5) >>> S = Sj.tospecdata() #Make spectrum ob diff --git a/pywafo/src/wafo/spectrum/models.py b/pywafo/src/wafo/spectrum/models.py index 9f90f22..9b08ac7 100644 --- a/pywafo/src/wafo/spectrum/models.py +++ b/pywafo/src/wafo/spectrum/models.py @@ -53,8 +53,6 @@ from numpy import (inf, atleast_1d, newaxis, any, minimum, maximum, array, #@Unr isfinite, mod, expm1, tanh, cosh, finfo, ones, ones_like, isnan, #@UnresolvedImport zeros_like, flatnonzero, sinc, hstack, vstack, real, flipud, clip) #@UnresolvedImport from dispersion_relation import w2k -#ppimport.enable() -#_wafospectrum = ppimport.ppimport('wafo.spectrum') from wafo.spectrum import SpecData1D sech = lambda x: 1.0 / cosh(x) @@ -638,8 +636,6 @@ class Jonswap(ModelSpectrum): def phi1(wi, h, g=9.81): ''' Factor transforming spectra to finite water depth spectra. - CALL: tr = phi1(w,h) - Input ----- w : arraylike @@ -743,7 +739,8 @@ class Tmaspec(Jonswap): phi1, Torsethaugen - References: + References + ---------- Buows, E., Gunther, H., Rosenthal, W., and Vincent, C.L. (1985) 'Similarity of the wind wave spectrum in finite depth water: 1 spectral form.' J. Geophys. Res., Vol 90, No. C1, pp 975-986 @@ -1361,8 +1358,8 @@ class Spreading(object): (Hasselman: spa ~= spb) (cos-2s) [6.97 9.77 4.06 -2.3 0 1.05 3 ] (Banner : spa ~= spb) (sech2) [2.61 2.28 1.3 -1.3 0.56 0.95 1.6] - Examples : - + Examples + -------- >>> import pylab as plb >>> D = Spreading('cos2s',s_a=10.0) @@ -1385,9 +1382,12 @@ class Spreading(object): >>> plb.close('all') - See also mkdspec, plotspec, spec2spec + See also + -------- + mkdspec, plotspec, spec2spec References + --------- Krogstad, H.E. and Barstow, S.F. (1999) "Directional Distributions in Ocean Wave Spectra" Proceedings of the 9th ISOPE Conference, Vol III, pp. 79-86 diff --git a/pywafo/src/wafo/transform/core.py b/pywafo/src/wafo/transform/core.py index 8ab3de5..d08b9ef 100644 --- a/pywafo/src/wafo/transform/core.py +++ b/pywafo/src/wafo/transform/core.py @@ -159,7 +159,7 @@ class TrData(WafoData, TrCommon): >>> g.sigma 5 >>> g.dat2gauss(1,2,3) - + [array([ 0.]), array([ 0.4]), array([ 0.6])] Check that the departure from a Gaussian model is zero >>> g.dist2gauss() < 1e-16 diff --git a/pywafo/src/wafo/transform/models.py b/pywafo/src/wafo/transform/models.py index c81c9db..7384abd 100644 --- a/pywafo/src/wafo/transform/models.py +++ b/pywafo/src/wafo/transform/models.py @@ -3,6 +3,7 @@ Transform Gaussian models ------------------------- TrHermite TrOchi +TrLinear ''' #------------------------------------------------------------------------------- # Name: transform.models @@ -84,7 +85,10 @@ class TrHermite(TrCommon): See also -------- - spec2skew, ochitr, lc2tr, dat2tr + SpecData1d.stats_nl + wafo.transform.TrOchi + wafo.objects.LevelCrossings.trdata + wafo.objects.TimeSeries.trdata References ---------- @@ -306,6 +310,11 @@ class TrLinear(TrCommon): See also -------- + TrOchi + TrHermite + SpecData1D.stats_nl + LevelCrossings.trdata + TimeSeries.trdata spec2skew, ochitr, lc2tr, dat2tr """ diff --git a/pywafo/src/wafo/transform/test/__init__.py b/pywafo/src/wafo/transform/test/__init__.py new file mode 100644 index 0000000..e69de29 diff --git a/pywafo/src/wafo/transform/test/test_models.py b/pywafo/src/wafo/transform/test/test_models.py new file mode 100644 index 0000000..ce3d66c --- /dev/null +++ b/pywafo/src/wafo/transform/test/test_models.py @@ -0,0 +1,47 @@ +from wafo.transform.models import TrHermite, TrOchi, TrLinear + +def test_trhermite(): + ''' + >>> std = 7./4 + >>> g = TrHermite(sigma=std, ysigma=std) + >>> g.dist2gauss() + 3.9858776379926808 + + >>> g.mean + 0.0 + >>> g.sigma + 1.75 + >>> g.dat2gauss([0,1,2,3]) + array([ 0.04654321, 1.03176393, 1.98871279, 2.91930895]) + + ''' +def test_trochi(): + ''' + >>> std = 7./4 + >>> g = TrOchi(sigma=std, ysigma=std) + >>> g.dist2gauss() + 5.9322684525265501 + >>> g.mean + 0.0 + >>> g.sigma + 1.75 + >>> g.dat2gauss([0,1,2,3]) + array([ 6.21927960e-04, 9.90237621e-01, 1.96075606e+00, + 2.91254576e+00]) + ''' +def test_trlinear(): + ''' + >>> std = 7./4 + >>> g = TrLinear(sigma=std, ysigma=std) + >>> g.dist2gauss() + 0.0 + >>> g.mean + 0.0 + >>> g.sigma + 1.75 + >>> g.dat2gauss([0,1,2,3]) + array([ 0., 1., 2., 3.]) + ''' +if __name__=='__main__': + import doctest + doctest.testmod() \ No newline at end of file diff --git a/pywafo/src/wafo/transform/test/test_trdata.py b/pywafo/src/wafo/transform/test/test_trdata.py new file mode 100644 index 0000000..b37e053 --- /dev/null +++ b/pywafo/src/wafo/transform/test/test_trdata.py @@ -0,0 +1,32 @@ +from wafo.transform import TrData + +def test_trdata(): + ''' + Construct a linear transformation model + >>> import numpy as np + >>> sigma = 5; mean = 1 + >>> u = np.linspace(-5,5); x = sigma*u+mean; y = u + >>> g = TrData(y,x) + >>> g.mean + array([ 1.]) + >>> g.sigma + array([ 5.]) + + >>> g = TrData(y,x,mean=1,sigma=5) + >>> g.mean + 1 + >>> g.sigma + 5 + >>> g.dat2gauss(1,2,3) + [array([ 0.]), array([ 0.4]), array([ 0.6])] + >>> g.dat2gauss([0,1,2,3]) + array([-0.2, 0. , 0.2, 0.4]) + + Check that the departure from a Gaussian model is zero + >>> g.dist2gauss() < 1e-16 + True + ''' + +if __name__=='__main__': + import doctest + doctest.testmod() \ No newline at end of file