diff --git a/.gitignore b/.gitignore
index 741aff5c..53f5bcc2 100644
--- a/.gitignore
+++ b/.gitignore
@@ -107,3 +107,6 @@ ENV/
 .idea/
 
 pip-wheel-metadata/
+
+# macOS DS_Store
+.DS_Store
diff --git a/examples/plot_oneway.py b/examples/plot_oneway.py
index 40ca6e47..06d8b68b 100644
--- a/examples/plot_oneway.py
+++ b/examples/plot_oneway.py
@@ -13,7 +13,8 @@ real dataset.
 
 import skfda
 from skfda.inference.anova import oneway_anova
-from skfda.representation import FDataGrid
+from skfda.representation import FDataGrid, FDataBasis
+from skfda.representation.basis import Fourier
 
 ################################################################################
 # *One-way ANOVA* (analysis of variance) is a test that can be used to
@@ -30,7 +31,7 @@ from skfda.representation import FDataGrid
 # of hips and knees from 39 different boys in a 20 point movement cycle.
 dataset = skfda.datasets.fetch_gait()
 fd_hip = dataset['data'].coordinates[0]
-fd_knee = dataset['data'].coordinates[1]
+fd_knee = dataset['data'].coordinates[1].to_basis(Fourier(n_basis=10))
 
 ################################################################################
 # Let's start with the first feature, the angle of the hip. The sample
@@ -49,8 +50,7 @@ fd_hip3 = fd_hip[26:39]
 fd_hip.plot(group=[0 if i < 13 else 1 if i < 26 else 39 for i in range(39)])
 
 means = [fd_hip1.mean(), fd_hip2.mean(), fd_hip3.mean()]
-fd_means = fd_hip.copy(data_matrix=[mean.data_matrix[0] for mean in means],
-                       dataset_label='Hip angle (means)')
+fd_means = skfda.concatenate(means)
 fig = fd_means.plot()
 
 ###############################################################################
@@ -74,7 +74,7 @@ print('p-value: ', p_val)
 
 ################################################################################
 # This was the simplest way to call this function. Let's see another example,
-# this time using knee angles.
+# this time using knee angles, this time with data in basis representation.
 fig = fd_knee.plot()
 
 ################################################################################
@@ -86,8 +86,7 @@ fd_knee3 = fd_knee[26:39]
 fd_knee.plot(group=[0 if i < 13 else 1 if i < 26 else 39 for i in range(39)])
 
 means = [fd_knee1.mean(), fd_knee2.mean(), fd_knee3.mean()]
-fd_means = fd_knee.copy(data_matrix=[mean.data_matrix[0] for mean in means],
-                        dataset_label='Knee angle (means)')
+fd_means = skfda.concatenate(means)
 fig = fd_means.plot()
 
 ################################################################################
@@ -104,7 +103,7 @@ fig = fd_means.plot()
 # sampling distribution of the statistic which is compared with the first
 # return to get the *p-value*.
 
-v_n, p_val, dist = oneway_anova(fd_knee1, fd_knee2, fd_knee3, n_reps=1500, p=2,
+v_n, p_val, dist = oneway_anova(fd_knee1, fd_knee2, fd_knee3, n_reps=1500,
                                 return_dist=True)
 
 print('Statistic: ', v_n)
diff --git a/examples/plot_oneway_synthetic.py b/examples/plot_oneway_synthetic.py
index 4c5583c6..a8844488 100644
--- a/examples/plot_oneway_synthetic.py
+++ b/examples/plot_oneway_synthetic.py
@@ -11,9 +11,13 @@ synthetic data.
 
 # sphinx_gallery_thumbnail_number = 2
 
+import numpy as np
+
 import skfda
-from skfda.inference.anova import oneway_anova
 from skfda.representation import FDataGrid
+from skfda.inference.anova import oneway_anova
+from skfda.datasets import make_gaussian_process
+from skfda.misc.covariances import WhiteNoise
 
 ################################################################################
 # *One-way ANOVA* (analysis of variance) is a test that can be used to
@@ -32,12 +36,6 @@ from skfda.representation import FDataGrid
 # test is to illustrate the differences in the results of the ANOVA method
 # when the covariance function of the brownian processes changes.
 
-import numpy as np
-
-from skfda.representation import FDataGrid
-from skfda.inference.anova import oneway_anova
-from skfda.datasets import make_gaussian_process
-
 ################################################################################
 # First, the means for the future processes are drawn.
 
@@ -69,8 +67,8 @@ groups[20:] = 'Sample 3'
 # differences between the means of each group should be clear, and the
 # p-value for the test should be near to zero.
 
-sigma = 0.01
-cov = np.identity(n_features) * sigma
+sigma2 = 0.01
+cov = WhiteNoise(variance=sigma2)
 
 fd1 = make_gaussian_process(n_samples, mean=m1, cov=cov,
                             n_features=n_features, random_state=1, start=start,
@@ -89,11 +87,12 @@ print("p-value: {:.3f}".format(p_val))
 # In the plot below we can see the simulated trajectories for each mean,
 # and the averages for each group.
 
-fd = fd1.concatenate(fd2.concatenate(fd3.concatenate()))
-fd.dataset_label = "Sample with $\sigma$ = {}, p-value = {:.3f}".format(
-    sigma, p_val)
-fd.plot(group=groups, legend=True, alpha=0.6)
-fd1.mean().concatenate(fd2.mean().concatenate(fd3.mean()).concatenate()).plot()
+fd = skfda.concatenate([fd1, fd2, fd3])
+fd_total = skfda.concatenate([fd.mean() for fd in [fd1, fd2,
+                                                               fd3]])
+fd_total.dataset_label = "Sample with $\sigma^2$ = {}, p-value = {:.3f}".format(
+    sigma2, p_val)
+fd_total.plot()
 
 ################################################################################
 # In the following, the same process will be followed incrementing sigma
@@ -103,9 +102,8 @@ fd1.mean().concatenate(fd2.mean().concatenate(fd3.mean()).concatenate()).plot()
 
 ################################################################################
 # Plot for :math:`\sigma = 1`:
-
-sigma = 1
-cov = np.identity(n_features) * sigma
+sigma2 = 0.1
+cov = WhiteNoise(variance=sigma2)
 
 fd1 = make_gaussian_process(n_samples, mean=m1, cov=cov,
                             n_features=n_features, random_state=1, start=t[0],
@@ -119,17 +117,18 @@ fd3 = make_gaussian_process(n_samples, mean=m3, cov=cov,
 
 _, p_val = oneway_anova(fd1, fd2, fd3, random_state=4)
 
-fd = fd1.concatenate(fd2.concatenate(fd3.concatenate()))
-fd.dataset_label = "Sample with $\sigma$ = {}, p-value = {:.3f}".format(
-    sigma, p_val)
-fd.plot(group=groups, legend=True, alpha=0.6)
-fd1.mean().concatenate(fd2.mean().concatenate(fd3.mean()).concatenate()).plot()
+fd = skfda.concatenate([fd1, fd2, fd3])
+fd_total = skfda.concatenate([fd.mean() for fd in [fd1, fd2,
+                                                               fd3]])
+fd_total.dataset_label = "Sample with $\sigma^2$ = {}, p-value = {:.3f}".format(
+    sigma2, p_val)
+fd_total.plot()
 
 ################################################################################
 # Plot for :math:`\sigma = 10`:
 
-sigma = 10
-cov = np.identity(n_features) * sigma
+sigma2 = 1
+cov = WhiteNoise(variance=sigma2)
 
 fd1 = make_gaussian_process(n_samples, mean=m1, cov=cov,
                             n_features=n_features, random_state=1, start=t[0],
@@ -143,11 +142,12 @@ fd3 = make_gaussian_process(n_samples, mean=m3, cov=cov,
 
 _, p_val = oneway_anova(fd1, fd2, fd3, random_state=4)
 
-fd = fd1.concatenate(fd2.concatenate(fd3.concatenate()))
-fd.dataset_label = "Sample with $\sigma$ = {}, p-value = {:.3f}".format(
-    sigma, p_val)
-fd.plot(group=groups, legend=True, alpha=0.6)
-fd1.mean().concatenate(fd2.mean().concatenate(fd3.mean()).concatenate()).plot()
+fd = skfda.concatenate([fd1, fd2, fd3])
+fd_total = skfda.concatenate([fd.mean() for fd in [fd1, fd2,
+                                                               fd3]])
+fd_total.dataset_label = "Sample with $\sigma^2$ = {}, p-value = {:.3f}".format(
+    sigma2, p_val)
+fd_total.plot()
 
 ################################################################################
 # **References:**
diff --git a/skfda/__init__.py b/skfda/__init__.py
index f1354bf9..c66f69d3 100644
--- a/skfda/__init__.py
+++ b/skfda/__init__.py
@@ -32,6 +32,7 @@ import errno as _errno
 from .representation import FData
 from .representation import FDataBasis
 from .representation import FDataGrid
+from .representation._functional_data import concatenate
 
 from . import representation, datasets, preprocessing, exploratory, misc, ml
 
diff --git a/skfda/inference/anova/anova_oneway.py b/skfda/inference/anova/anova_oneway.py
index 03ab8558..8ebc1137 100644
--- a/skfda/inference/anova/anova_oneway.py
+++ b/skfda/inference/anova/anova_oneway.py
@@ -1,37 +1,38 @@
 import numpy as np
 from sklearn.utils import check_random_state
 
-from skfda.misc.metrics import norm_lp
-from skfda.representation import FData, FDataGrid, FDataBasis
+from skfda import concatenate
+from skfda.misc.metrics import lp_distance
+from skfda.representation import FData, FDataGrid
 from skfda.datasets import make_gaussian_process
 
 
 def v_sample_stat(fd, weights, p=2):
     r"""
     Calculates a statistic that measures the variability between groups of
-    samples in a :class:`skfda.representation.grid.FDataGrid` object.
+    samples in a :class:`skfda.representation.FData` object.
 
     The statistic defined as below is calculated between all the samples in a
-    :class:`skfda.representation.grid.FDataGrid` object with a given set of
-    weights, and the desired :math:`L_p` norm.
+    :class:`skfda.representation.FData` object with a given set of
+    weights.
 
-    Let :math:`\{f_i\}_{i=1}^k` be a set of samples in a FDataGrid object.
+    Let :math:`\{f_i\}_{i=1}^k` be a set of samples in a FData object.
     Let :math:`\{w_j\}_{j=1}^k` be a set of weights, where :math:`w_i` is
     related to the sample :math:`f_i` for :math:`i=1,\dots,k`.
     The statistic is defined as:
 
     .. math::
-        V_n = \sum_{i<j}^kw_i\|f_i-f_j\|^p
+        V_n = \sum_{i<j}^kw_i\|f_i-f_j\|^2
 
     This statistic is defined in Cuevas[1].
 
     Args:
-         fd (FDataGrid): Object containing all the samples for which we want
+         fd (FData): Object containing all the samples for which we want
             to calculate the statistic.
          weights (list of int): Weights related to each sample. Each
             weight is expected to appear in the same position as its
-            corresponding sample in the FDataGrid object.
-         p (int, optional): p of the lp norm. Must be greater or equal
+            corresponding sample in the FData object.
+        p (int, optional): p of the lp norm. Must be greater or equal
             than 1. If p='inf' or p=np.inf it is used the L infinity metric.
             Defaults to 2.
 
@@ -68,47 +69,43 @@ def v_sample_stat(fd, weights, p=2):
         Analysis*, 47:111-112, 02 2004
     """
 
+    weights = np.asarray(weights)
     if not isinstance(fd, FData):
         raise ValueError("Argument type must inherit FData.")
     if len(weights) != fd.n_samples:
         raise ValueError("Number of weights must match number of samples.")
-    if isinstance(fd, FDataBasis):
-        raise NotImplementedError("Not implemented for FDataBasis objects.")
 
-    n = fd.n_samples
-    v_n = 0
-    for j in range(n):
-        for i in range(j):
-            v_n += weights[i] * norm_lp(fd[i] - fd[j], p=p) ** p
-    return v_n
+    t_ind = np.tril_indices(fd.n_samples, -1)
+    coef = weights[t_ind[1]]
+    return np.sum(coef * lp_distance(fd[t_ind[0]], fd[t_ind[1]], p=p) ** p)
 
 
 def v_asymptotic_stat(fd, weights, p=2):
     r"""
     Calculates a statistic that measures the variability between groups of
-    samples in a :class:`skfda.representation.grid.FDataGrid` object.
+    samples in a :class:`skfda.representation.FData` object.
 
     The statistic defined as below is calculated between all the samples in a
-    :class:`skfda.representation.grid.FDataGrid` object with a given set of
-    weights, and the desired :math:`L_p` norm.
+    :class:`skfda.representation.FData` object with a given set of
+    weights.
 
-    Let :math:`\{f_i\}_{i=1}^k` be a set of samples in a FDataGrid object.
+    Let :math:`\{f_i\}_{i=1}^k` be a set of samples in a FData object.
     Let :math:`\{w_j\}_{j=1}^k` be a set of weights, where :math:`w_i` is
     related to the sample :math:`f_i` for :math:`i=1,\dots,k`.
     The statistic is defined as:
 
     .. math::
-        \sum_{i<j}^k\|f_i-f_j\sqrt{\cfrac{w_i}{w_j}}\|^p
+        \sum_{i<j}^k\|f_i-f_j\sqrt{\cfrac{w_i}{w_j}}\|^2
 
     This statistic is defined in Cuevas[1].
 
     Args:
-         fd (FDataGrid): Object containing all the samples for which we want
+         fd (FData): Object containing all the samples for which we want
             to calculate the statistic.
          weights (list of int): Weights related to each sample. Each
             weight is expected to appear in the same position as its
-            corresponding sample in the FDataGrid object.
-         p (int, optional): p of the lp norm. Must be greater or equal
+            corresponding sample in the FData object.
+        p (int, optional): p of the lp norm. Must be greater or equal
             than 1. If p='inf' or p=np.inf it is used the L infinity metric.
             Defaults to 2.
 
@@ -116,7 +113,7 @@ def v_asymptotic_stat(fd, weights, p=2):
         The value of the statistic.
 
     Raises:
-        ValueError.
+        ValueError
 
     Examples:
 
@@ -144,36 +141,35 @@ def v_asymptotic_stat(fd, weights, p=2):
         anova test for functional data". *Computational Statistics  Data
         Analysis*, 47:111-112, 02 2004
     """
+    weights = np.asarray(weights)
     if not isinstance(fd, FData):
         raise ValueError("Argument type must inherit FData.")
     if len(weights) != fd.n_samples:
         raise ValueError("Number of weights must match number of samples.")
-    if isinstance(fd, FDataBasis):
-        raise NotImplementedError("Not implemented for FDataBasis objects.")
+    if np.count_nonzero(weights) != len(weights):
+        raise ValueError("All weights must be non-zero.")
 
-    n = fd.n_samples
-    v = 0
-    for j in range(n):
-        for i in range(j):
-            v += norm_lp(
-                fd[i] - fd[j] * np.sqrt(weights[i] / weights[j]), p=p) ** p
-    return v
+    t_ind = np.tril_indices(fd.n_samples, -1)
+    coef = np.sqrt(weights[t_ind[1]] / weights[t_ind[0]])
+    left_fd = fd[t_ind[1]]
+    if isinstance(fd, FDataGrid):
+        right_fd = coef[:, None, np.newaxis] * fd[t_ind[0]]
+    else:
+        right_fd = fd[t_ind[0]].times(coef)
+    return np.sum(lp_distance(left_fd, right_fd, p=p) ** p)
 
 
-def _anova_bootstrap(fd_grouped, n_reps, p=2, random_state=None):
+def _anova_bootstrap(fd_grouped, n_reps, random_state=None, p=2):
 
     n_groups = len(fd_grouped)
-    assert n_groups > 0
+    if n_groups < 2:
+        raise ValueError("At least two groups must be passed in fd_grouped.")
 
-    # Creating list with all the sample points
-    list_sample = [fd.sample_points[0].tolist() for fd in fd_grouped]
-    # Checking that the all the entries in the list are the same
-    if not list_sample.count(list_sample[0]) == len(list_sample):
-        raise ValueError("All FDataGrid passed must have the same sample "
-                         "points.")
+    for fd in fd_grouped[1:]:
+        if not np.array_equal(fd.domain_range, fd_grouped[0].domain_range):
+            raise ValueError("Domain range must match for every FData in "
+                             "fd_grouped.")
 
-    sample_points = fd_grouped[0].sample_points
-    m = len(sample_points[0])  # Number of points in the grid
     start, stop = fd_grouped[0].domain_range[0]
 
     sizes = [fd.n_samples for fd in fd_grouped]  # List with sizes of each group
@@ -181,12 +177,17 @@ def _anova_bootstrap(fd_grouped, n_reps, p=2, random_state=None):
     # Estimating covariances for each group
     k_est = [fd.cov().data_matrix[0, ..., 0] for fd in fd_grouped]
 
+    # Number of sample points for gaussian processes have to match the features
+    # of the covariances.
+    n_features = k_est[0].shape[0]
+
     # Instance a random state object in case random_state is an int
     random_state = check_random_state(random_state)
 
     # Simulating n_reps observations for each of the n_groups gaussian processes
-    sim = [make_gaussian_process(n_reps, n_features=m, start=start, stop=stop,
-                                 cov=k_est[i], random_state=random_state)
+    sim = [make_gaussian_process(n_reps, n_features=n_features, start=start,
+                                 stop=stop, cov=k_est[i],
+                                 random_state=random_state)
            for i in range(n_groups)]
     v_samples = np.empty(n_reps)
     for i in range(n_reps):
@@ -195,7 +196,7 @@ def _anova_bootstrap(fd_grouped, n_reps, p=2, random_state=None):
     return v_samples
 
 
-def oneway_anova(*args, n_reps=2000, p=2, return_dist=False, random_state=None):
+def oneway_anova(*args, n_reps=2000, return_dist=False, random_state=None, p=2):
     r"""
     Performs one-way functional ANOVA.
 
@@ -233,16 +234,15 @@ def oneway_anova(*args, n_reps=2000, p=2, return_dist=False, random_state=None):
             procedure. Defaults to 2000 (This value may change in future
             versions).
 
-        p (int, optional): p of the lp norm. Must be greater or equal
-            than 1. If p='inf' or p=np.inf it is used the L infinity metric.
-            Defaults to 2.
-
         return_dist (bool, optional): Flag to indicate if the function should
-        return a
-            numpy.array with the sampling distribution simulated.
+            return a numpy.array with the sampling distribution simulated.
 
         random_state (optional): Random state.
 
+        p (int, optional): p of the lp norm. Must be greater or equal
+            than 1. If p='inf' or p=np.inf it is used the L infinity metric.
+            Defaults to 2.
+
     Returns:
         Value of the sample statistic, p-value and sampling distribution of
         the simulated asymptotic statistic.
@@ -263,8 +263,6 @@ def oneway_anova(*args, n_reps=2000, p=2, return_dist=False, random_state=None):
         >>> fd1, fd2, fd3 = fd[:13], fd[13:26], fd[26:]
         >>> oneway_anova(fd1, fd2, fd3, random_state=RandomState(42))
         (179.52499999999998, 0.602)
-        >>> oneway_anova(fd1, fd2, fd3, p=1, random_state=RandomState(42))
-        (67.27499999999999, 0.0)
         >>> _, _, dist = oneway_anova(fd1, fd2, fd3, n_reps=3,
         ...     random_state=RandomState(42),
         ...     return_dist=True)
@@ -284,25 +282,38 @@ def oneway_anova(*args, n_reps=2000, p=2, return_dist=False, random_state=None):
         raise ValueError("Argument type must inherit FData.")
     if n_reps < 1:
         raise ValueError("Number of simulations must be positive.")
-    if any(isinstance(fd, FDataBasis) for fd in args):
-        raise NotImplementedError("Not implemented for FDataBasis objects.")
 
     fd_groups = args
-    # Creating list with all the sample points
-    list_sample = [fd.sample_points[0].tolist() for fd in fd_groups]
-    # Checking that the all the entries in the list are the same
-    if not list_sample.count(list_sample[0]) == len(list_sample):
-        raise ValueError("All FDataGrid passed must have the same sample "
-                         "points.")
-
-    fd_means = fd_groups[0].mean()
-    for fd in fd_groups[1:]:
-        fd_means = fd_means.concatenate(fd.mean())
+    if not all([isinstance(fd, type(fd_groups[0])) for fd in fd_groups[1:]]):
+        raise TypeError('Found mixed FData types in arguments.')
 
+    for fd in fd_groups[1:]:
+        if not np.array_equal(fd.domain_range, fd_groups[0].domain_range):
+            raise ValueError("Domain range must match for every FData passed.")
+
+    if isinstance(fd_groups[0], FDataGrid):
+        # Creating list with all the sample points
+        list_sample = [fd.sample_points[0].tolist() for fd in fd_groups]
+        # Checking that the all the entries in the list are the same
+        if not list_sample.count(list_sample[0]) == len(list_sample):
+            raise ValueError("All FDataGrid passed must have the same sample "
+                             "points.")
+    else:  # If type is FDataBasis, check same basis
+        list_basis = [fd.basis for fd in fd_groups]
+        if not list_basis.count(list_basis[0]) == len(list_basis):
+            raise NotImplementedError("Not implemented for FDataBasis with "
+                                      "different basis.")
+
+    # FData where each sample is the mean of each group
+    fd_means = concatenate([fd.mean() for fd in fd_groups])
+
+    # Base statistic
     vn = v_sample_stat(fd_means, [fd.n_samples for fd in fd_groups], p=p)
 
-    simulation = _anova_bootstrap(fd_groups, n_reps, p=p,
-                                  random_state=random_state)
+    # Computing sampling distribution
+    simulation = _anova_bootstrap(fd_groups, n_reps,
+                                  random_state=random_state, p=p)
+
     p_value = np.sum(simulation > vn) / len(simulation)
 
     if return_dist:
diff --git a/skfda/misc/covariances.py b/skfda/misc/covariances.py
index 4c4d154f..806f6897 100644
--- a/skfda/misc/covariances.py
+++ b/skfda/misc/covariances.py
@@ -401,10 +401,29 @@ class Exponential(Covariance):
         y = _transform_to_2d(y)
 
         x_y = _squared_norms(x, y)
-
         return self.variance * np.exp(-np.sqrt(x_y) / (self.length_scale))
 
     def to_sklearn(self):
         """Convert it to a sklearn kernel, if there is one"""
         return (self.variance *
                 sklearn_kern.Matern(length_scale=self.length_scale, nu=0.5))
+
+
+class WhiteNoise(Covariance):
+    """Gaussian covariance function."""
+
+    _latex_formula = (r"K(x,y)= \left\{ \begin{array}{lc} \sigma^2, & x = y "
+                      r"\\ 0, & x \neq y\\ \end{array} \right.")
+
+    _parameters = [("variance", r"\sigma^2")]
+
+    def __init__(self, *, variance: float = 1.):
+        self.variance = variance
+
+    def __call__(self, x, y):
+        x = _transform_to_2d(x)
+        return self.variance * np.eye(x.shape[0])
+
+    def to_sklearn(self):
+        """Convert it to a sklearn kernel, if there is one"""
+        return sklearn_kern.WhiteKernel(noise_level=self.variance)
diff --git a/skfda/misc/metrics.py b/skfda/misc/metrics.py
index e05a99d7..ba70aa5e 100644
--- a/skfda/misc/metrics.py
+++ b/skfda/misc/metrics.py
@@ -6,7 +6,7 @@ from ..preprocessing.registration import normalize_warping, ElasticRegistration
 from ..preprocessing.registration._warping import _normalize_scale
 from ..preprocessing.registration.elastic import SRSF
 from ..representation import FData
-from ..representation import FDataGrid
+from ..representation import FDataGrid, FDataBasis
 
 
 def _cast_to_grid(fdata1, fdata2, eval_points=None, _check=True, **kwargs):
@@ -209,7 +209,7 @@ def pairwise_distance(distance, **kwargs):
     return pairwise
 
 
-def norm_lp(fdatagrid, p=2, p2=2):
+def norm_lp(fdata, p=2, p2=2):
     r"""Calculate the norm of all the samples in a FDataGrid object.
 
     For each sample sample f the Lp norm is defined as:
@@ -244,7 +244,7 @@ def norm_lp(fdatagrid, p=2, p2=2):
 
 
     Args:
-        fdatagrid (FDataGrid): FDataGrid object.
+        fdata (FData): FData object.
         p (int, optional): p of the lp norm. Must be greater or equal
             than 1. If p='inf' or p=np.inf it is used the L infinity metric.
             Defaults to 2.
@@ -279,30 +279,40 @@ def norm_lp(fdatagrid, p=2, p2=2):
     if not (p == 'inf' or np.isinf(p)) and p < 1:
         raise ValueError(f"p must be equal or greater than 1.")
 
-    if fdatagrid.dim_codomain > 1:
-        if p2 == 'inf':
-            p2 = np.inf
-        data_matrix = np.linalg.norm(fdatagrid.data_matrix, ord=p2, axis=-1,
-                                     keepdims=True)
-    else:
-        data_matrix = np.abs(fdatagrid.data_matrix)
+    if isinstance(fdata, FDataBasis):
+        if fdata.dim_codomain > 1 or p != 2:
+            raise NotImplementedError
 
-    if p == 'inf' or np.isinf(p):
+        start, end = fdata.domain_range[0]
+        integral = scipy.integrate.quad_vec(
+            lambda x: np.power(np.abs(fdata(x)), p), start, end)
+        res = np.sqrt(integral[0]).flatten()
 
-        if fdatagrid.dim_domain == 1:
-            res = np.max(data_matrix[..., 0], axis=1)
+    else:
+        if fdata.dim_codomain > 1:
+            if p2 == 'inf':
+                p2 = np.inf
+            data_matrix = np.linalg.norm(fdata.data_matrix, ord=p2, axis=-1,
+                                         keepdims=True)
         else:
-            res = np.array([np.max(sample) for sample in data_matrix])
+            data_matrix = np.abs(fdata.data_matrix)
 
-    elif fdatagrid.dim_domain == 1:
+        if p == 'inf' or np.isinf(p):
 
-        # Computes the norm, approximating the integral with Simpson's rule.
-        res = scipy.integrate.simps(data_matrix[..., 0] ** p,
-                                    x=fdatagrid.sample_points) ** (1 / p)
+            if fdata.dim_domain == 1:
+                res = np.max(data_matrix[..., 0], axis=1)
+            else:
+                res = np.array([np.max(sample) for sample in data_matrix])
 
-    else:
-        # Needed to perform surface integration
-        return NotImplemented
+        elif fdata.dim_domain == 1:
+
+            # Computes the norm, approximating the integral with Simpson's rule.
+            res = scipy.integrate.simps(data_matrix[..., 0] ** p,
+                                        x=fdata.sample_points) ** (1 / p)
+
+        else:
+            # Needed to perform surface integration
+            return NotImplemented
 
     if len(res) == 1:
         return res[0]
diff --git a/skfda/representation/_functional_data.py b/skfda/representation/_functional_data.py
index 5a1a0294..9ba2c101 100644
--- a/skfda/representation/_functional_data.py
+++ b/skfda/representation/_functional_data.py
@@ -1009,3 +1009,34 @@ class FData(ABC, pandas.api.extensions.ExtensionArray):
         first, *others = to_concat
 
         return first.concatenate(*others)
+
+
+def concatenate(objects, as_coordinates=False):
+    """
+    Join samples from an iterable of similar FData objects.
+
+    Joins samples of FData objects if they have the same
+    dimensions and sampling points.
+    Args:
+        objects (list of :obj:`FDataBasis`): Objects to be concatenated.
+        as_coordinates (boolean, optional):  If False concatenates as
+                new samples, else, concatenates the other functions as
+                new components of the image. Defaults to False.
+    Returns:
+        :obj:`FData`: FData object with the samples from the
+        original objects.
+    Raises:
+        ValueError: In case the provided list of FData objects is
+        empty.
+    Todo:
+        By the moment, only unidimensional objects are supported in basis
+        representation.
+    """
+    objects = iter(objects)
+    first = next(objects, None)
+
+    if not first:
+        raise ValueError("At least one FData object must be provided "
+                         "to concatenate.")
+
+    return first.concatenate(*objects, as_coordinates=as_coordinates)

Concatenate function for FData objects of the same type
In this moment to concatenate `FData` objects it is necessary to call the `concatenate()` method. For example, let `fd_list`be a list of `FDataGrid` objects. To get the concatenation of the items in the list is necessary to execute:
```python
fd_c = fd_list[0].concatenate(*fd_list[1:])
```

It would be interesting to implement a function that receives a list of `FData` objects of the same type and returns the concatenated `FData` object. For example:
```python
fd_c = fdata_concatenate(fd_list)
```
diff --git a/tests/test_basis.py b/tests/test_basis.py
index ca6b29b3..26ef8199 100644
--- a/tests/test_basis.py
+++ b/tests/test_basis.py
@@ -1,5 +1,7 @@
 from skfda.representation.basis import (Basis, FDataBasis, Constant, Monomial,
                                         BSpline, Fourier)
+from skfda.representation.grid import FDataGrid
+from skfda import concatenate
 import unittest
 
 import numpy as np
@@ -432,6 +434,19 @@ class TestBasis(unittest.TestCase):
                                         [-120, -18, -60],
                                         [-48, 0, 48]])
 
+    def test_concatenate(self):
+        sample1 = np.arange(0, 10)
+        sample2 = np.arange(10, 20)
+        fd1 = FDataGrid([sample1]).to_basis(Fourier(n_basis=5))
+        fd2 = FDataGrid([sample2]).to_basis(Fourier(n_basis=5))
+
+        fd = concatenate([fd1, fd2])
+
+        np.testing.assert_equal(fd.n_samples, 2)
+        np.testing.assert_equal(fd.dim_codomain, 1)
+        np.testing.assert_equal(fd.dim_domain, 1)
+        np.testing.assert_array_equal(fd.coefficients, np.concatenate(
+            [fd1.coefficients, fd2.coefficients]))
 
 if __name__ == '__main__':
     print()
diff --git a/tests/test_covariances.py b/tests/test_covariances.py
index 8eccd806..a4e29024 100644
--- a/tests/test_covariances.py
+++ b/tests/test_covariances.py
@@ -13,9 +13,8 @@ class TestsSklearn(unittest.TestCase):
 
     def _test_compare_sklearn(self, cov: skfda.misc.covariances.Covariance):
         cov_sklearn = cov.to_sklearn()
-
         cov_matrix = cov(self.x, self.x)
-        cov_sklearn_matrix = cov_sklearn(self.x, self.x)
+        cov_sklearn_matrix = cov_sklearn(self.x)
 
         np.testing.assert_array_almost_equal(cov_matrix, cov_sklearn_matrix)
 
@@ -62,3 +61,10 @@ class TestsSklearn(unittest.TestCase):
                     cov = skfda.misc.covariances.Exponential(
                         variance=variance, length_scale=length_scale)
                     self._test_compare_sklearn(cov)
+
+    def test_white_noise(self):
+
+        for variance in [1, 2]:
+            with self.subTest(variance=variance):
+                cov = skfda.misc.covariances.WhiteNoise(variance=variance)
+                self._test_compare_sklearn(cov)
diff --git a/tests/test_grid.py b/tests/test_grid.py
index 2026cefa..4213b145 100644
--- a/tests/test_grid.py
+++ b/tests/test_grid.py
@@ -3,7 +3,7 @@ import unittest
 import scipy.stats.mstats
 
 import numpy as np
-from skfda import FDataGrid
+from skfda import FDataGrid, concatenate
 from skfda.exploratory import stats
 
 
@@ -99,6 +99,22 @@ class TestFDataGrid(unittest.TestCase):
         fd = fd1.concatenate(fd2, as_coordinates=True)
         np.testing.assert_equal(None, fd.axes_labels)
 
+    def test_concatenate(self):
+        sample1 = np.arange(0, 10)
+        sample2 = np.arange(10, 20)
+        fd1 = FDataGrid([sample1])
+        fd2 = FDataGrid([sample2])
+
+        fd1.axes_labels = ["x", "y"]
+        fd = concatenate([fd1, fd2])
+
+        np.testing.assert_equal(fd.n_samples, 2)
+        np.testing.assert_equal(fd.dim_codomain, 1)
+        np.testing.assert_equal(fd.dim_domain, 1)
+        np.testing.assert_array_equal(fd.data_matrix[..., 0], [sample1,
+                                                               sample2])
+        np.testing.assert_array_equal(fd1.axes_labels, fd.axes_labels)
+
     def test_coordinates(self):
         fd1 = FDataGrid([[1, 2, 3, 4], [2, 3, 4, 5]])
         fd1.axes_labels = ["x", "y"]
diff --git a/tests/test_oneway_anova.py b/tests/test_oneway_anova.py
index 4e34d329..f7c2ed87 100644
--- a/tests/test_oneway_anova.py
+++ b/tests/test_oneway_anova.py
@@ -1,7 +1,9 @@
 import unittest
 import numpy as np
+import pytest
 
 from skfda.representation import FDataGrid
+from skfda.representation.basis import Fourier
 from skfda.datasets import fetch_gait
 from skfda.inference.anova import oneway_anova, v_asymptotic_stat, \
     v_sample_stat
@@ -26,6 +28,9 @@ class OnewayAnovaTests(unittest.TestCase):
             v_asymptotic_stat(1, [1])
         with self.assertRaises(ValueError):
             v_asymptotic_stat(FDataGrid([0]), [0, 1])
+        with self.assertRaises(ValueError):
+            v_asymptotic_stat(FDataGrid([[1, 1, 1], [1, 1, 1]]), [0, 0])
+
 
     def test_v_stats(self):
         n_features = 50
@@ -36,28 +41,38 @@ class OnewayAnovaTests(unittest.TestCase):
         m3 = [3 for _ in range(n_features)]
         fd = FDataGrid([m1, m2, m3], sample_points=t)
         self.assertEqual(v_sample_stat(fd, weights), 7.0)
-        self.assertEqual(v_sample_stat(fd, weights, p=1), 5.0)
+        self.assertAlmostEqual(v_sample_stat(fd.to_basis(Fourier(n_basis=5)),
+                                             weights), 7.0)
         res = (1 - 2 * np.sqrt(1 / 2)) ** 2 + (1 - 3 * np.sqrt(1 / 3)) ** 2 \
                                             + (2 - 3 * np.sqrt(2 / 3)) ** 2
         self.assertAlmostEqual(v_asymptotic_stat(fd, weights), res)
-        res = abs(1 - 2 * np.sqrt(1 / 2)) + abs(1 - 3 * np.sqrt(1 / 3))\
-                                          + abs(2 - 3 * np.sqrt(2 / 3))
-        self.assertAlmostEqual(v_asymptotic_stat(fd, weights, p=1), res)
+        self.assertAlmostEqual(v_asymptotic_stat(fd.to_basis(Fourier(
+            n_basis=5)), weights), res)
 
     def test_asymptotic_behaviour(self):
         dataset = fetch_gait()
         fd = dataset['data'].coordinates[1]
-        fd1 = fd[0:13]
-        fd2 = fd[13:26]
-        fd3 = fd[26:39]
+        fd1 = fd[0:5]
+        fd2 = fd[5:10]
+        fd3 = fd[10:15]
+
+        n_little_sim = 10
+
+        sims = np.array([oneway_anova(fd1, fd2, fd3, n_reps=500)[1] for _ in
+                         range(n_little_sim)])
+        little_sim = np.mean(sims)
+        big_sim = oneway_anova(fd1, fd2, fd3, n_reps=2000)[1]
+        self.assertAlmostEqual(little_sim, big_sim, delta=0.05)
 
-        n_little_sim = 50
+        fd = fd.to_basis(Fourier(n_basis=5))
+        fd1 = fd[0:5]
+        fd2 = fd[5:10]
 
-        sims = np.array([oneway_anova(fd1, fd2, fd3, n_reps=2000)[1] for _ in
+        sims = np.array([oneway_anova(fd1, fd2, n_reps=500)[1] for _ in
                          range(n_little_sim)])
         little_sim = np.mean(sims)
-        big_sim = oneway_anova(fd1, fd2, fd3, n_reps=50000)[1]
-        self.assertAlmostEqual(little_sim, big_sim, delta=0.01)
+        big_sim = oneway_anova(fd1, fd2, n_reps=2000)[1]
+        self.assertAlmostEqual(little_sim, big_sim, delta=0.05)
 
 
 if __name__ == '__main__':

pip install -e .
parse_log_pytest
requirements.txt
pytest pytest-cov
apt-get update apt-get install -y gcc pip3 install --upgrade pip cython numpy
3.6
requirements.txt
pytest --no-header -rA --tb=line --color=no -p no:cacheprovider -W ignore::DeprecationWarning
attrs==22.2.0
certifi==2021.5.30
coverage==6.2
cycler==0.11.0
Cython==3.0.12
findiff==0.10.0
importlib-metadata==4.8.3
iniconfig==1.1.1
joblib==1.1.1
kiwisolver==1.3.1
matplotlib==3.3.4
mpldatacursor==0.7.1
mpmath==1.3.0
multimethod==1.5
numpy==1.19.5
packaging==21.3
pandas==1.1.5
Pillow==8.4.0
pluggy==1.0.0
py==1.11.0
pyparsing==3.1.4
pytest==7.0.1
pytest-cov==4.0.0
python-dateutil==2.9.0.post0
pytz==2025.2
rdata==0.3
scikit-datasets==0.2.0
-e git+https://github.com/GAA-UAM/scikit-fda.git@6727106eca0f51d1fd8142495e1c722adf50f2dc#egg=scikit_fda
scikit-learn==0.24.2
scipy==1.5.4
six==1.17.0
sklearn==0.0
sympy==1.9
threadpoolctl==3.1.0
tomli==1.2.3
typing_extensions==4.1.1
xarray==0.16.2
zipp==3.6.0

name: scikit-fda
channels:
  - defaults
  - https://repo.anaconda.com/pkgs/main
  - https://repo.anaconda.com/pkgs/r
  - conda-forge
dependencies:
  - _libgcc_mutex=0.1=main
  - _openmp_mutex=5.1=1_gnu
  - ca-certificates=2025.2.25=h06a4308_0
  - certifi=2021.5.30=py36h06a4308_0
  - ld_impl_linux-64=2.40=h12ee557_0
  - libffi=3.3=he6710b0_2
  - libgcc-ng=11.2.0=h1234567_1
  - libgomp=11.2.0=h1234567_1
  - libstdcxx-ng=11.2.0=h1234567_1
  - ncurses=6.4=h6a678d5_0
  - openssl=1.1.1w=h7f8727e_0
  - python=3.6.13=h12debd9_1
  - readline=8.2=h5eee18b_0
  - setuptools=58.0.4=py36h06a4308_0
  - sqlite=3.45.3=h5eee18b_0
  - tk=8.6.14=h39e8969_0
  - wheel=0.37.1=pyhd3eb1b0_0
  - xz=5.6.4=h5eee18b_1
  - zlib=1.2.13=h5eee18b_1
  - pip:
      - attrs==22.2.0
      - coverage==6.2
      - cycler==0.11.0
      - cython==3.0.12
      - findiff==0.10.0
      - importlib-metadata==4.8.3
      - iniconfig==1.1.1
      - joblib==1.1.1
      - kiwisolver==1.3.1
      - matplotlib==3.3.4
      - mpldatacursor==0.7.1
      - mpmath==1.3.0
      - multimethod==1.5
      - numpy==1.19.5
      - packaging==21.3
      - pandas==1.1.5
      - pillow==8.4.0
      - pip==21.3.1
      - pluggy==1.0.0
      - py==1.11.0
      - pyparsing==3.1.4
      - pytest==7.0.1
      - pytest-cov==4.0.0
      - python-dateutil==2.9.0.post0
      - pytz==2025.2
      - rdata==0.3
      - scikit-datasets==0.2.0
      - scikit-learn==0.24.2
      - scipy==1.5.4
      - six==1.17.0
      - sklearn==0.0
      - sympy==1.9
      - threadpoolctl==3.1.0
      - tomli==1.2.3
      - typing-extensions==4.1.1
      - xarray==0.16.2
      - zipp==3.6.0
prefix: /opt/conda/envs/scikit-fda

tests/test_basis.py::TestBasis::test_basis_basis_inprod tests/test_basis.py::TestBasis::test_basis_bspline_product tests/test_basis.py::TestBasis::test_basis_constant_product tests/test_basis.py::TestBasis::test_basis_fdatabasis_inprod tests/test_basis.py::TestBasis::test_basis_fourier_product tests/test_basis.py::TestBasis::test_basis_gram_matrix tests/test_basis.py::TestBasis::test_basis_inner_matrix tests/test_basis.py::TestBasis::test_basis_monomial_product tests/test_basis.py::TestBasis::test_basis_product_generic tests/test_basis.py::TestBasis::test_comutativity_inprod tests/test_basis.py::TestBasis::test_concatenate tests/test_basis.py::TestBasis::test_fdatabasis__add__ tests/test_basis.py::TestBasis::test_fdatabasis__mul__ tests/test_basis.py::TestBasis::test_fdatabasis__mul__2 tests/test_basis.py::TestBasis::test_fdatabasis__sub__ tests/test_basis.py::TestBasis::test_fdatabasis_derivative_bspline tests/test_basis.py::TestBasis::test_fdatabasis_derivative_constant tests/test_basis.py::TestBasis::test_fdatabasis_derivative_fourier tests/test_basis.py::TestBasis::test_fdatabasis_derivative_monomial tests/test_basis.py::TestBasis::test_fdatabasis_fdatabasis_inprod tests/test_basis.py::TestBasis::test_fdatabasis_times_fdatabasis_fdatabasis tests/test_basis.py::TestBasis::test_fdatabasis_times_fdatabasis_int tests/test_basis.py::TestBasis::test_fdatabasis_times_fdatabasis_list tests/test_basis.py::TestBasis::test_from_data_cholesky tests/test_basis.py::TestBasis::test_from_data_qr tests/test_covariances.py::TestsSklearn::test_exponential tests/test_covariances.py::TestsSklearn::test_gaussian tests/test_covariances.py::TestsSklearn::test_linear tests/test_covariances.py::TestsSklearn::test_polynomial tests/test_covariances.py::TestsSklearn::test_white_noise tests/test_grid.py::TestFDataGrid::test_add tests/test_grid.py::TestFDataGrid::test_concatenate tests/test_grid.py::TestFDataGrid::test_concatenate_coordinates tests/test_grid.py::TestFDataGrid::test_coordinates tests/test_grid.py::TestFDataGrid::test_copy_equals tests/test_grid.py::TestFDataGrid::test_gmean tests/test_grid.py::TestFDataGrid::test_init tests/test_grid.py::TestFDataGrid::test_mean tests/test_grid.py::TestFDataGrid::test_slice tests/test_oneway_anova.py::OnewayAnovaTests::test_oneway_anova_args tests/test_oneway_anova.py::OnewayAnovaTests::test_v_stats tests/test_oneway_anova.py::OnewayAnovaTests::test_v_stats_args
tests/test_oneway_anova.py::OnewayAnovaTests::test_asymptotic_behaviour
skfda/misc/covariances.py skfda/misc/metrics.py examples/plot_oneway_synthetic.py skfda/inference/anova/anova_oneway.py .gitignore examples/plot_oneway.py skfda/representation/_functional_data.py skfda/__init__.py
skfda/misc/covariances.py skfda/misc/metrics.py examples/plot_oneway_synthetic.py skfda/inference/anova/anova_oneway.py .gitignore examples/plot_oneway.py skfda/representation/_functional_data.py skfda/__init__.py