def heavsineFunction(dim=1024, normalize=True, heavs=0):
    """
    Generate the Heavsine function in a range of 0 to 1, with dim points.

    Parameters
    ----------
    dim: int
        Dimension of the signal.

    normalize: bool, optional
        It is True by default. This parameter normalize the data values 
        in a range of 0 to 1 with a function present in 
        ``statsWaveletFilt.miscellaneous``.

    heavs: int or float, optional
        It is 0 by default. This parameter, called * heavs * is the number of 
        discontinuities in the heavens characteristic signal shown by 
        Donoho [1] with 0 the signal will be the original, used in [1].

    Returns
    -------
    tuple:
        [0] 1-D array-like, coordinates in X axis and [1] 1-D array-like,
        coordinates in Y axis

    References
    ----------
    .. [1] DONOHO, D. L.; JOHNSTONE, I. M. Ideal spatial adaptation via
           wavelet shrinkage. Biometrika, v. 81, p. 425–455, 1994.
    """
    import numpy as np
    import statsWaveletFilt.miscellaneous as misc
    linspace = np.linspace
    sin = np.sin
    signal = np.sign
    pi = np.pi
    x = linspace(0, 1, dim)
    if heavs == 0:
        y = 4 * sin(4 * pi * x) - signal(x - 0.3) - signal(0.72 - x)
    else:
        y = 4 * sin(4 * pi * x)
        for i in range(heavs):
            if i % 2:
                y -= signal(x - np.random.random())
            else:
                y += signal(x - np.random.random())
    if normalize:
        y_nor = misc.normalizeData(y)
    else:
        y_nor = y
    return x, y_nor
