A method for performing multiple access in wireless OFDM cellular systems considering both space and frequency domains, comprising at least one base station equipped with a large number of antennas according to a two-dimensional rectangular array and a number M of cell users, said rectangular array comprising N 1 antenna elements along one axis with a regular spacing d x and N 2 antenna elements along a perpendicular axis with a regular spacing d y, said cell users being characterized by angles ( θ,φ) in a spherical coordinate system, characterized in that in order to achieve orthogonal multiple access the method comprises using an orthogonal space-frequency processing unit for: - selecting a grid spacing ( Δu,Δv) in the (u, v) domain through the equations: Δ ⁢ u = λ N 1 d x, Δ ⁢ v = λ N 2 d y where λ denotes the wavelength of the system operating frequency, and the (u, v) domain is obtained from the following transformation of the spherical angles ( θ,φ): u = sin θ cos φ; v = sin θ sin φ - discretizing the (u, v) domain according to the following expressions: u k = k ⋅ Δ ⁢ u, k = 0,1,..., N 1 − 1 v l = l ⋅ Δ ⁢ v, l = 0,1,..., N 2 − 1, where the visible region of the grid corresponds to those values enclosed within a unit circle according to the equation: u 2 + v 2 ≤ 1; - constructing a set of signals S T [ k,l,f] according to the following expression: S T k l f = { S i f; ∀ k l = k i l i, i = 0,1,..., M − 1 0; otherwise, where S i [ f] denote the frequency contents of the complex baseband signal aimed at user i, and ( k i,l i) represent the estimated user locations in the discretized (u, v) grid; - calculating the time-domain excitations A T [ n,m,t] for the antenna elements in the array given by coordinates ( nd x,md y), for generation of the downlink transmit signals according to the following expression: A T n m t = 1 N c N 1 N 2 ∑ f = 0 N c − 1 ∑ k = 0 N 1 − 1 ∑ l = 0 N 2 − 1 S T k l f exp j 2 ⁢ π N c ft exp − j 2 ⁢ π N 1 nk exp − j 2 ⁢ π N 2 ml, where N c denotes the number of subcarriers in the system bandwidth; and - obtaining the frequency contents S R [ k,l,f] of the complex baseband signals received from the M users in the uplink, by applying the following transformation over a set of received signals A R [ n,m,t] corresponding to each of the antenna elements: S R k l f = ∑ t = 0 N c − 1 ∑ n = 0 N 1 − 1 ∑ m = 0 N 2 − 1 A R n m t exp − j 2 ⁢ π N c ft exp j 2 ⁢ π N 1 nk exp j 2 ⁢ π N 2 ml. - selecting a grid spacing ( Δu,Δv) in the (u, v) domain through the equations: Δ ⁢ u = λ N 1 d x, Δ ⁢ v = λ N 2 d y where λ denotes the wavelength of the system operating frequency, and the (u, v) domain is obtained from the following transformation of the spherical angles ( θ,φ): u = sin θ cos φ; v = sin θ sin φ - discretizing the (u, v) domain according to the following expressions: u k = k ⋅ Δ ⁢ u, k = 0,1,..., N 1 − 1 v l = l ⋅ Δ ⁢ v, l = 0,1,..., N 2 − 1, where the visible region of the grid corresponds to those values enclosed within a unit circle according to the equation: u 2 + v 2 ≤ 1; - constructing a set of signals S T [ k,l,f] according to the following expression: S T k l f = { S i f; ∀ k l = k i l i, i = 0,1,..., M − 1 0; otherwise, where S i [ f] denote the frequency contents of the complex baseband signal aimed at user i, and ( k i,l i) represent the estimated user locations in the discretized (u, v) grid; - calculating the time-domain excitations A T [ n,m,t] for the antenna elements in the array given by coordinates ( nd x,md y), for generation of the downlink transmit signals according to the following expression: A T n m t = 1 N c N 1 N 2 ∑ f = 0 N c − 1 ∑ k = 0 N 1 − 1 ∑ l = 0 N 2 − 1 S T k l f exp j 2 ⁢ π N c ft exp − j 2 ⁢ π N 1 nk exp − j 2 ⁢ π N 2 ml, where N c denotes the number of subcarriers in the system bandwidth; and - obtaining the frequency contents S R [ k,l,f] of the complex baseband signals received from the M users in the uplink, by applying the following transformation over a set of received signals A R [ n,m,t] corresponding to each of the antenna elements: S R k l f = ∑ t = 0 N c − 1 ∑ n = 0 N 1 − 1 ∑ m = 0 N 2 − 1 A R n m t exp − j 2 ⁢ π N c ft exp j 2 ⁢ π N 1 nk exp j 2 ⁢ π N 2 ml.