diff --git a/pywafo/src/wafo/fig.py b/pywafo/src/wafo/fig.py
index 77bcbb7..5dacb32 100644
--- a/pywafo/src/wafo/fig.py
+++ b/pywafo/src/wafo/fig.py
@@ -41,10 +41,10 @@ import win32con
 import msvcrt
 import numpy
 
-__all__ = ['close','cycle','hide','keep','maximize','minimize','pile','restore','stack','tile']
+__all__ = ['close', 'cycle', 'hide', 'keep', 'maximize', 'minimize', 'pile', 'restore', 'stack', 'tile']
 
 # Figure format strings to recognize in window title
-_FIG_FORMATS = ('Figure','TVTK Scene','Chaco Plot Window: Figure')
+_FIG_FORMATS = ('Figure', 'TVTK Scene', 'Chaco Plot Window: Figure')
 _SCREENSIZE = None
 
 def _getScreenSize(wnds):
@@ -98,7 +98,7 @@ def _findTopWindows(wantedTitle=None):
     if wantedTitle == None:
         return topWindows
     else:
-        return [(hwnd,windowTxt) for hwnd, windowTxt in topWindows if windowTxt.startswith(wantedTitle)]
+        return [(hwnd, windowTxt) for hwnd, windowTxt in topWindows if windowTxt.startswith(wantedTitle)]
 
 def findallfigs():
     '''
@@ -120,10 +120,10 @@ def findallfigs():
 def _figparse(*args):
     figs = []
     for arg in args:
-        if  isinstance(arg,(list,tuple,set)):
+        if  isinstance(arg, (list, tuple, set)):
             for val in arg:
                 figs.append(int(val))
-        elif isinstance(arg,int):
+        elif isinstance(arg, int):
             figs.append(arg)
         elif arg == 'all':
             figs = 'all'
@@ -132,7 +132,7 @@ def _figparse(*args):
             raise TypeError('Only integers arguments accepted!')
             #raise TypeError('Unrecognized argument type (%s)!'%type(arg))
 
-    if len(figs)==0 or figs =='all':
+    if len(figs) == 0 or figs == 'all':
         figs = findallfigs()
     return figs
 
@@ -145,12 +145,12 @@ def _fig2wnd(figs):
     for fig in figs:
         for format_ in _FIG_FORMATS:
             winTitle = format_ + ' %d' % fig
-            hwnd = FindWindow(None,winTitle)
+            hwnd = FindWindow(None, winTitle)
             if not hwnd == 0:
                 wnd_handles.append(hwnd)
     return wnd_handles
 
-def _show_figure(figs,cmdshow):
+def _show_figure(figs, cmdshow):
     ''' sets the specified figure's show state.
 
     @param figs: vector for figure numbers
@@ -186,13 +186,13 @@ def _show_figure(figs,cmdshow):
     for fig in figs:
         for format_ in _FIG_FORMATS:
             winTitle = format_ + ' %d' % fig
-            hwnd = FindWindow(None,winTitle)
+            hwnd = FindWindow(None, winTitle)
             if not hwnd == 0:
                 #ShowWindow(hwnd,cmdshow)
                 BringWindowToTop(hwnd)
-                ShowWindow(hwnd,cmdshow)
+                ShowWindow(hwnd, cmdshow)
 
-def _show_windows(wnds,cmdshow):
+def _show_windows(wnds, cmdshow):
     ''' sets the specified window's show state.
 
     @param wnds: list of window handles numbers
@@ -228,7 +228,7 @@ def _show_windows(wnds,cmdshow):
         if not hwnd == 0:
             #ShowWindow(hwnd,cmdshow)
             BringWindowToTop(hwnd)
-            ShowWindow(hwnd,cmdshow)
+            ShowWindow(hwnd, cmdshow)
 
 def keep(*figs):
     ''' Keeps figure windows of your choice and closes the rest.
@@ -255,15 +255,15 @@ def keep(*figs):
     '''
     figs2keep = []
     for fig in figs:
-        if  isinstance(fig,(list,tuple,set)):
+        if  isinstance(fig, (list, tuple, set)):
             for val in fig:
                 figs2keep.append(int(val))
-        elif isinstance(fig,int):
+        elif isinstance(fig, int):
             figs2keep.append(fig)
         else:
             raise TypeError('Only integers arguments accepted!')
 
-    if len(figs2keep)>0:
+    if len(figs2keep) > 0:
         allfigs = set(findallfigs())
 
 # Remove figure handles in the "keep" list
@@ -476,16 +476,16 @@ def pile(*figs):
     figlist = _figparse(*figs)
     wnds = _fig2wnd(figlist)
     numfigs = len(wnds)
-    if numfigs>0:
+    if numfigs > 0:
         pos = _getScreenSize(wnds)
-        pos[3] = int(pos[3]/2)
-        pos[2] = int(pos[2]/2.5)
-        pos[1] = int(pos[3]/2)
-        pos[0] = int(pos[2]/2)
+        pos[3] = int(pos[3] / 2)
+        pos[2] = int(pos[2] / 2.5)
+        pos[1] = int(pos[3] / 2)
+        pos[0] = int(pos[2] / 2)
         BringWindowToTop = win32gui.BringWindowToTop
         MoveWindow = win32gui.MoveWindow
         for wnd in wnds:
-            MoveWindow(wnd, pos[0], pos[1], pos[2], pos[3],1)
+            MoveWindow(wnd, pos[0], pos[1], pos[2], pos[3], 1)
             BringWindowToTop(wnd)
 
 
@@ -521,17 +521,17 @@ def stack(*figs):
     figlist = _figparse(*figs)
     wnds = _fig2wnd(figlist)
     numfigs = len(wnds)
-    if numfigs>0:
+    if numfigs > 0:
         screenpos = _getScreenSize(wnds)
 
-        maxfigs = numpy.fix(screenpos[3]/20)
+        maxfigs = numpy.fix(screenpos[3] / 20)
 
-        if (numfigs>maxfigs):            # figure limit check
-            print(' More than %d requested '% maxfigs)
+        if (numfigs > maxfigs):            # figure limit check
+            print(' More than %d requested ' % maxfigs)
             return
         BringWindowToTop = win32gui.BringWindowToTop
         MoveWindow = win32gui.MoveWindow
-        GetWindowRect =  win32gui.GetWindowRect
+        GetWindowRect = win32gui.GetWindowRect
 #
 #
 #% tile figures by postiion
@@ -543,13 +543,13 @@ def stack(*figs):
             pos[3] -= pos[1]
             pos[2] -= pos[0]
             #print('[x, y, w, h] = ', pos)
-            ypos = screenpos[1]+iy*20 #int(screenpos[3] - iy*20 -pos[3] -70) # figure location (row)
-            xpos = int(iy*5 + 15+screenpos[0])     # figure location (column)
-            MoveWindow(wnds[iy],xpos,ypos, pos[2],pos[3],1)
+            ypos = screenpos[1] + iy * 20 #int(screenpos[3] - iy*20 -pos[3] -70) # figure location (row)
+            xpos = int(iy * 5 + 15 + screenpos[0])     # figure location (column)
+            MoveWindow(wnds[iy], xpos, ypos, pos[2], pos[3], 1)
             BringWindowToTop(wnds[iy])
 
 
-def tile(*figs,**kwds):
+def tile(*figs, **kwds):
     ''' Tile figure windows.
 
     Parameters
@@ -587,17 +587,17 @@ def tile(*figs,**kwds):
 
     nfigs = len(wnds);        # Number of windows.
 
-    if nfigs>0:
-        nfigspertile = kwds.get('pairs',nfigs)
+    if nfigs > 0:
+        nfigspertile = kwds.get('pairs', nfigs)
 
         ceil = numpy.ceil
         sqrt = numpy.sqrt
         maximum = numpy.maximum
 
-        nlayers = int(ceil(nfigs/nfigspertile))
+        nlayers = int(ceil(nfigs / nfigspertile))
 
         nh = int(ceil(sqrt(nfigspertile)))   # Number of figures horisontally.
-        nv = int(ceil(nfigspertile/nh));     # Number of figures vertically.
+        nv = int(ceil(nfigspertile / nh));     # Number of figures vertically.
 
 
         nh = maximum(nh, 2);
@@ -622,15 +622,15 @@ def tile(*figs,**kwds):
 #% 4 - Window vertical size
 #%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
         
-        hspc   = 10            # Horisontal space.
+        hspc = 10            # Horisontal space.
         topspc = 20;            # Space above top figure.
         medspc = 10;            # Space between figures.
         botspc = 20;            # Space below bottom figure.
 
         #print('scrwid = %d' % scrwid)
-        figwid = ( scrwid - (nh+1)*hspc )/nh
+        figwid = (scrwid - (nh + 1) * hspc) / nh
         #print('figwid = %d' % figwid)
-        fighgt = ( scrhgt - (topspc+botspc) - (nv-1)*medspc )/nv;
+        fighgt = (scrhgt - (topspc + botspc) - (nv - 1) * medspc) / nv;
 
         figwid = int(numpy.round(figwid))
         fighgt = int(numpy.round(fighgt))
@@ -642,19 +642,19 @@ def tile(*figs,**kwds):
         for unused_ix in xrange(nlayers):
             for row in xrange(nv):
                 for col in xrange(nh):
-                    if  (row)*nh + col < nfigspertile :      
-                        if idx<nfigs:
-                            figlft = int(screenpos[0]+(col+1)*hspc + col*figwid)
-                            figtop = int(screenpos[1]+topspc + row*(fighgt + medspc))
+                    if  (row) * nh + col < nfigspertile :      
+                        if idx < nfigs:
+                            figlft = int(screenpos[0] + (col + 1) * hspc + col * figwid)
+                            figtop = int(screenpos[1] + topspc + row * (fighgt + medspc))
                             #figpos = [ figlft figtop figwid fighgt ];    # Figure position.
                             #fighnd = FindWindow(0,'Figure %d' % figs[idx]) # Figure handle.
                             fighnd = wnds[idx]
-                            MoveWindow(fighnd, figlft,figtop, figwid, fighgt,1);     # Set position.
+                            MoveWindow(fighnd, figlft, figtop, figwid, fighgt, 1);     # Set position.
                             BringWindowToTop(fighnd)
                             #figure(figs[idx]);                              # Raise figure.
                         idx += 1
 
-def cycle(*figs,**kwds):
+def cycle(*figs, **kwds):
     ''' Cycle through figure windows.
 
     Parameters
@@ -697,9 +697,9 @@ def cycle(*figs,**kwds):
     wnds = _fig2wnd(figlist)
 
     numfigs = len(wnds)
-    if numfigs>0:
-        maximize = kwds.get('maximize',False)
-        pairs = kwds.get('pairs',1)
+    if numfigs > 0:
+        maximize = kwds.get('maximize', False)
+        pairs = kwds.get('pairs', 1)
 
         if maximize or pairs == None:
             nfigspercycle = 1
@@ -713,18 +713,18 @@ def cycle(*figs,**kwds):
         i = 0;
         escape_key = chr(27)
         backspace_key = chr(8)
-        while 0<=i and i<numfigs:
+        while 0 <= i and i < numfigs:
 
             if maximize:
                 cmdshow = win32con.SW_SHOWMAXIMIZED
             else:
                 cmdshow = win32con.SW_SHOWNORMAL
 
-            iu = min(i+nfigspercycle,numfigs)
+            iu = min(i + nfigspercycle, numfigs)
             wnd = wnds[i:iu]
             _show_windows(wnd, cmdshow)
 
-            if i + nfigspercycle-1< numfigs:
+            if i + nfigspercycle - 1 < numfigs:
                 print('Press escape to quit, backspace to display previous figure(s) and any other key to display next figure(s)');
 
             #time.sleep(0.5)
@@ -735,9 +735,9 @@ def cycle(*figs,**kwds):
                 _show_windows(wnd, win32con.SW_RESTORE)
 
 
-            if B==backspace_key: # Back space
+            if B == backspace_key: # Back space
                 i -= nfigspercycle
-            elif B==escape_key:
+            elif B == escape_key:
                 break
             else:
                 i += nfigspercycle
@@ -747,9 +747,7 @@ def cycle(*figs,**kwds):
         _show_windows(wnds, win32con.SW_SHOWNORMAL)
 
 
-
-
-if __name__=='__main__':
+if __name__ == '__main__':
     import doctest
     doctest.testmod()
 
@@ -810,4 +808,4 @@ if __name__=='__main__':
 #    ShowWindow = prototype3(("ShowWindow", windll.user32),paramflags)
 #    import win32gui
 #    h = win32gui.FindWindow(None,'PyLab')
-#    win32gui.ShowWindow(h,)
\ No newline at end of file
+#    win32gui.ShowWindow(h,)
diff --git a/pywafo/src/wafo/kdetools.py b/pywafo/src/wafo/kdetools.py
index 6a5a864..a07a780 100644
--- a/pywafo/src/wafo/kdetools.py
+++ b/pywafo/src/wafo/kdetools.py
@@ -19,7 +19,7 @@ import scipy.optimize as optimize
 from wafo.misc import meshgrid
 from wafo.wafodata import WafoData
 
-from dctpack import dct
+from dctpack import dct, dctn, idctn
 
 import copy
 import numpy as np
@@ -63,7 +63,220 @@ def sphere_volume(d, r=1.0):
     Chapman and Hall, pp 105
     """
     return (r ** d) * 2.0 * pi ** (d / 2.0) / (d * gamma(d / 2.0))
+class KDEgauss(object):
+    """ Kernel-Density Estimator base class.
+
+    Parameters
+    ----------
+    data : (# of dims, # of data)-array
+        datapoints to estimate from
+    hs : array-like (optional) 
+        smooting parameter vector/matrix.
+        (default compute from data using kernel.get_smoothing function)
+    alpha : real scalar (optional)
+        sensitivity parameter               (default 0 regular KDE)
+        A good choice might be alpha = 0.5 ( or 1/D)
+        alpha = 0      Regular  KDE (hs is constant)
+        0 < alpha <= 1 Adaptive KDE (Make hs change)  
+
+
+    Members
+    -------
+    d : int
+        number of dimensions
+    n : int
+        number of datapoints
+
+    Methods
+    -------
+    kde.eval_grid_fast(x0, x1,..., xd) : array
+        evaluate the estimated pdf on meshgrid(x0, x1,..., xd)
+    kde(x0, x1,..., xd) : array
+        same as kde.eval_grid_fast(x0, x1,..., xd)
+    """
+
+    def __init__(self, data, hs=None, kernel=None, alpha=0.0, xmin=None, xmax=None, inc=128):
+        self.dataset = atleast_2d(data)
+        self.hs = hs
+        self.kernel = kernel if kernel else Kernel('gauss')
+        self.alpha = alpha
+        self.xmin = xmin
+        self.xmax = xmax
+        self.inc = inc
+        self.initialize()
+
+    def initialize(self):
+        self.d, self.n = self.dataset.shape
+        self._set_xlimits()
+        self._initialize()
+        
+    def _initialize(self):
+        self._compute_smoothing()
+        
+    def _compute_smoothing(self):
+        """Computes the smoothing matrix
+        """
+        get_smoothing = self.kernel.get_smoothing
+        h = self.hs
+        if h is None:
+            h = get_smoothing(self.dataset)
+        h = np.atleast_1d(h)
+        hsiz = h.shape
+    
+        if (len(hsiz) == 1) or (self.d == 1):
+            if max(hsiz) == 1:
+                h = h * np.ones(self.d)
+            else:
+                h.shape = (self.d,) # make sure it has the correct dimension
+          
+            # If h negative calculate automatic values
+            ind, = np.where(h <= 0)
+            for i in ind.tolist(): # 
+                h[i] = get_smoothing(self.dataset[i])
+            deth = h.prod()
+            self.inv_hs = linalg.diag(1.0 / h)
+        else: #fully general smoothing matrix
+            deth = linalg.det(h)
+            if deth <= 0:
+                raise ValueError('bandwidth matrix h must be positive definit!')
+            self.inv_hs = linalg.inv(h)
+        self.hs = h
+        self._norm_factor = deth * self.n
+        
+    def _set_xlimits(self):
+        amin = self.dataset.min(axis= -1)
+        amax = self.dataset.max(axis= -1)
+        iqr = iqrange(self.dataset, axis= -1)
+        sigma = np.minimum(np.std(self.dataset, axis= -1, ddof=1), iqr / 1.34)
+        #xyzrange = amax - amin
+        #offset = xyzrange / 4.0
+        offset = 2 * sigma
+        if self.xmin is None:
+            self.xmin = amin - offset
+        else:
+            self.xmin = self.xmin * np.ones((self.d,1))
+        if self.xmax is None:
+            self.xmax = amax + offset
+        else:
+            self.xmax = self.xmax * np.ones((self.d,1))
+            
+    def eval_grid_fast(self, *args, **kwds):
+        """Evaluate the estimated pdf on a grid.
+
+        Parameters
+        ----------
+        arg_0,arg_1,... arg_d-1 : vectors
+            Alternatively, if no vectors is passed in then
+             arg_i = linspace(self.xmin[i], self.xmax[i], self.inc)
+        output : string optional
+            'value' if value output
+            'data' if object output
+
+        Returns
+        -------
+        values : array-like
+            The values evaluated at meshgrid(*args).
+
+        """
+        if len(args) == 0:
+            args = []
+            for i in range(self.d):
+                args.append(np.linspace(self.xmin[i], self.xmax[i], self.inc))
+        self.args = args
+        return self._eval_grid_fun(self._eval_grid_fast, *args, **kwds)
+       
+    def _eval_grid_fast(self, *args, **kwds):
+        X = np.vstack(args)
+        d, inc = X.shape
+        dx = X[:, 1] - X[:, 0]
+        R = X.max(axis=-1)- X.min(axis=-1)
+        
+        t_star = (self.hs/R)**2
+        I = (np.asfarray(np.arange(0, inc))*pi)**2 
+        In = []
+        
+        for i in range(d):
+            In.append(I * t_star[i] * 0.5)
+        
+        Inc = meshgrid(*In) if d > 1 else In
+        
+        kw = np.zeros((inc,)*d)
+        for i in range(d):
+            kw += exp(-Inc[i])
+        y = kwds.get('y', 1.0)
+        d, n = self.dataset.shape
+        # Find the binned kernel weights, c.
+        c = gridcount(self.dataset, X, y=y)/n
+        # Perform the convolution.
+        at = dctn(c) * kw
+        z = idctn(at)*at.size/np.prod(R)
+        return z*(z>0.0)
+        #ix = (slice(0, inc),)*d
+        #return z[ix] * (z[ix] > 0.0)
+    
+    def _eval_grid_fun(self, eval_grd, *args, **kwds):
+        output = kwds.pop('output', 'value')
+        f = eval_grd(*args, **kwds)
+        if output == 'value':
+            return f
+        else:
+            titlestr = 'Kernel density estimate (%s)' % self.kernel.name
+            kwds2 = dict(title=titlestr)
+            kwds2['plot_kwds'] = dict(plotflag=1)
+            kwds2.update(**kwds)
+            args = self.args
+            if self.d == 1:
+                args = args[0]
+            wdata = WafoData(f, args, **kwds2)
+            if self.d > 1:
+                PL = np.r_[10:90:20, 95, 99, 99.9]
+                try:
+                    ql = qlevels(f, p=PL)
+                    wdata.clevels = ql
+                    wdata.plevels = PL
+                except:
+                    pass
+            return wdata 
 
+    def _check_shape(self, points):
+        points = atleast_2d(points)
+        d, m = points.shape
+        if d != self.d:
+            if d == 1 and m == self.d:
+                # points was passed in as a row vector
+                points = np.reshape(points, (self.d, 1))
+            else:
+                msg = "points have dimension %s, dataset has dimension %s" % (d,
+                    self.d)
+                raise ValueError(msg)
+        return points   
+    def eval_points(self, points, **kwds):
+        """Evaluate the estimated pdf on a set of points.
+
+        Parameters
+        ----------
+        points : (# of dimensions, # of points)-array
+            Alternatively, a (# of dimensions,) vector can be passed in and
+            treated as a single point.
+
+        Returns
+        -------
+        values : (# of points,)-array
+            The values at each point.
+
+        Raises
+        ------
+        ValueError if the dimensionality of the input points is different than
+        the dimensionality of the KDE.
+        """
+
+        points = self._check_shape(points)
+        return self._eval_points(points, **kwds)
+    
+    def _eval_points(self, points, **kwds):
+        pass
+
+    __call__ = eval_grid_fast
 class _KDE(object):
     """ Kernel-Density Estimator base class.
 
@@ -341,7 +554,7 @@ class TKDE(_KDE):
     def __init__(self, data, hs=None, kernel=None, alpha=0.0, xmin=None,
                  xmax=None, inc=128, L2=None):
         self.L2 = L2
-        _KDE.__init__(self, data, hs, kernel, alpha, xmin, xmax, inc)
+        super(TKDE, self).__init__(data, hs, kernel, alpha, xmin, xmax, inc)
     
     def _initialize(self):
         self._check_xmin()
@@ -574,7 +787,7 @@ class KDE(_KDE):
     """
 
     def __init__(self, data, hs=None, kernel=None, alpha=0.0, xmin=None, xmax=None, inc=128):
-        _KDE.__init__(self, data, hs, kernel, alpha, xmin, xmax, inc)
+        super(KDE, self).__init__(data, hs, kernel, alpha, xmin, xmax, inc)
             
     def _initialize(self):
         self._compute_smoothing()
@@ -638,11 +851,14 @@ class KDE(_KDE):
         Xn = np.dot(self.inv_hs, np.vstack(Xnc))
         
         # Obtain the kernel weights.        
+        kw = self.kernel(Xn) 
+        norm_fact = (kw.sum()*dx.prod()*self.n)
+        norm_fact2 = (self._norm_factor * self.kernel.norm_factor(d, self.n))
+        
+        kw = kw/norm_fact
         r = kwds.get('r', 0)
-        if r==0:
-            kw = self.kernel(Xn) / (self._norm_factor * self.kernel.norm_factor(d, self.n))
-        else:
-            kw = np.vstack(Xnc) ** r * self.kernel(Xn) / (self._norm_factor * self.kernel.norm_factor(d, self.n))
+        if r!=0:
+            kw *= np.vstack(Xnc) ** r 
         kw.shape = shape0
         kw = np.fft.ifftshift(kw)
         fftn = np.fft.fftn
@@ -1434,7 +1650,8 @@ class Kernel(object):
             
             c = gridcount(A[dim], xa)
             N = len(set(A[dim]))
-            a = dct(c/c.sum(), norm=None)
+            #a = dct(c/c.sum(), norm=None)
+            a = dct(c/len(A[dim]), norm=None)
             
             #% now compute the optimal bandwidth^2 using the referenced method
             I = np.asfarray(np.arange(1, inc))**2 
@@ -2443,7 +2660,7 @@ def gridcount(data, X, y=1):
     >>> import numpy as np
     >>> import wafo.kdetools as wk
     >>> import pylab as plb
-    >>> N     = 20;
+    >>> N = 20;
     >>> data  = np.random.rayleigh(1,N)
     >>> x = np.linspace(0,max(data)+1,50)  
     >>> dx = x[1]-x[0]  
@@ -2551,8 +2768,8 @@ def kde_demo1():
     import scipy.stats as st
     x = np.linspace(-4, 4, 101)
     x0 = x / 2.0
-    data = np.random.normal(loc=0, scale=1.0, size=7) #rndnorm(0,1,7,1);
-    kernel = Kernel('gaus')
+    data = np.random.normal(loc=0, scale=1.0, size=7) 
+    kernel = Kernel('gauss')
     hs = kernel.hns(data)
     hVec = [hs / 2, hs, 2 * hs]
        
@@ -2742,7 +2959,59 @@ def kreg_demo1(hs=None, fast=False, fun='hisj'):
     
     print(kreg.tkde.tkde.inv_hs)
     print(kreg.tkde.tkde.hs) 
-   
+
+def kde_gauss_demo(n=50000):
+    '''
+    KDEDEMO Demonstrate the KDEgauss
+
+    KDEDEMO1 shows the true density (dotted) compared to KDE based on 7
+    observations (solid) and their individual kernels (dashed) for 3
+    different values of the smoothing parameter, hs.
+    '''
+    
+    import scipy.stats as st
+    #x = np.linspace(-4, 4, 101)
+    #data = np.random.normal(loc=0, scale=1.0, size=n) 
+    #data = np.random.exponential(scale=1.0, size=n)
+#    n1 = 128
+#    I = (np.arange(n1)*pi)**2 *0.01*0.5
+#    kw = exp(-I)
+#    pylab.plot(idctn(kw))
+#    return
+    dist = st.norm
+    dist = st.expon
+    data = dist.rvs(loc=0, scale=1.0, size=n)
+    d, N = np.atleast_2d(data).shape
+    
+    if d==1:
+        plot_options = [dict(color='red'), dict(color='green'), dict(color='black')]
+    else:
+        plot_options = [dict(colors='red'), dict(colors='green'), dict(colors='black')]
+    
+    pylab.figure(1)
+    kde0 = KDE(data, kernel=Kernel('gauss', 'hisj'))
+    f0 = kde0.eval_grid_fast(output='plot', ylab='Density')
+    f0.plot(**plot_options[0])
+    
+    kde1 = TKDE(data, kernel=Kernel('gauss', 'hisj'), L2=0)
+    f1 = kde1.eval_grid_fast(output='plot', ylab='Density')
+    f1.plot(**plot_options[1])
+    
+    kde2 = KDEgauss(data)
+    f2 = kde2(output='plot', ylab='Density')
+    x = f2.args
+    f2.plot(**plot_options[2])
+    
+    fmax = dist.pdf(x, 0, 1).max() 
+    if d==1:
+        pylab.plot(x, dist.pdf(x, 0, 1), 'k:')
+        pylab.axis([x.min(), x.max(), 0, fmax])
+    pylab.show()
+    print(fmax/f2.data.max())
+    format_ = ''.join(('%g, ')*d)
+    format_ = 'hs0=%s hs2=%s' % (format_,format_)
+    print(format_ % tuple(kde0.hs.tolist()+kde2.hs.tolist()))
+        
 def test_docstrings():
     import doctest
     doctest.testmod()
@@ -2750,4 +3019,5 @@ def test_docstrings():
 if __name__ == '__main__':
     #test_docstrings()
     #kde_demo2()
-    kreg_demo1()
\ No newline at end of file
+    #kreg_demo1(fast=True)
+    kde_gauss_demo()
\ No newline at end of file
diff --git a/pywafo/src/wafo/wafodata.py b/pywafo/src/wafo/wafodata.py
index 062cec1..2440903 100644
--- a/pywafo/src/wafo/wafodata.py
+++ b/pywafo/src/wafo/wafodata.py
@@ -190,7 +190,7 @@ class Plotter_1d(object):
     def plot(self, wdata, *args, **kwds):
         plotflag = kwds.pop('plotflag', False)
         if plotflag:
-            h1 = self._plot(plotflag, wdata, **kwds)
+            h1 = self._plot(plotflag, wdata, *args, **kwds)
         else:
             if isinstance(wdata.args, (list, tuple)):
                 args1 = tuple((wdata.args)) + (wdata.data,) + args