A method for functional magnetic resonance imaging (fMRI) of a subject, the method comprising: applying with an MRI system, a balanced steady state free precession (b-SSFP) sequence to a target area in the subject; acquiring with the MRI system, image data of the target area in the subject using a randomly undersampled variable density spiral trajectory, wherein the VDS trajectory includes a random stack of VDS trajectories, wherein a change in a radial direction versus a change in a rate in an angular direction of the VDS trajectories is determined by a number of interleaves and an effective field of view, wherein the field of view follows an exponential function and an angle between each interleaf of the number of interleaves is randomly distributed; and producing an image of the target area in the subject, using compressed sensing reconstruction, based on the acquired image data, wherein an Archimedean spiral trajectory with constraints defined in the following equations is used to define the 3D random stack of VDS trajectories: Δk r Δθ = N 2 πFOV e FOV e = FOV 0 exp − k r ak rmax exp − k z bk zmax 2 2 where: the Archimedean spiral trajectory is described in polar coordinates as: k = k r e jθ where k = k x + ik y k r = k x 2 + k y 2 and θ = tan − 1 k y k x Δk r is the change rate in the radial direction Δθ is the change rate in the angular direction FOV e is the effective field of view in x and y FOV 0 is the desired field of view k rmax is the maximum value of k r which is one half of the resolution in x and y k z is the location for each k z - slice and k zmax is the number of k z -slices N is the number of interleaves parameters a and b control the variance of the exponential function along k r and k z direction and determine the acceleration factor, and wherein the number of interleaves on each k z - slice is designed to decrease with the increase of k z. applying with an MRI system, a balanced steady state free precession (b-SSFP) sequence to a target area in the subject; acquiring with the MRI system, image data of the target area in the subject using a randomly undersampled variable density spiral trajectory, wherein the VDS trajectory includes a random stack of VDS trajectories, wherein a change in a radial direction versus a change in a rate in an angular direction of the VDS trajectories is determined by a number of interleaves and an effective field of view, wherein the field of view follows an exponential function and an angle between each interleaf of the number of interleaves is randomly distributed; and producing an image of the target area in the subject, using compressed sensing reconstruction, based on the acquired image data, wherein an Archimedean spiral trajectory with constraints defined in the following equations is used to define the 3D random stack of VDS trajectories: Δk r Δθ = N 2 πFOV e FOV e = FOV 0 exp − k r ak rmax exp − k z bk zmax 2 2 where: the Archimedean spiral trajectory is described in polar coordinates as: k = k r e jθ where k = k x + ik y k r = k x 2 + k y 2 and θ = tan − 1 k y k x Δk r is the change rate in the radial direction Δθ is the change rate in the angular direction FOV e is the effective field of view in x and y FOV 0 is the desired field of view k rmax is the maximum value of k r which is one half of the resolution in x and y k z is the location for each k z - slice and k zmax is the number of k z -slices N is the number of interleaves parameters a and b control the variance of the exponential function along k r and k z direction and determine the acceleration factor, and wherein the number of interleaves on each k z - slice is designed to decrease with the increase of k z. the Archimedean spiral trajectory is described in polar coordinates as: k = k r e jθ where k = k x + ik y k r = k x 2 + k y 2 and θ = tan − 1 k y k x Δk r is the change rate in the radial direction Δθ is the change rate in the angular direction FOV e is the effective field of view in x and y FOV 0 is the desired field of view k rmax is the maximum value of k r which is one half of the resolution in x and y k z is the location for each k z - slice and k zmax is the number of k z -slices N is the number of interleaves parameters a and b control the variance of the exponential function along k r and k z direction and determine the acceleration factor, and wherein the number of interleaves on each k z - slice is designed to decrease with the increase of k z.