A widespread method for gyrokinetic simulations relies on the discretization of distribution functions with numerical markers. The electromagnetic field is typically represented on a grid by finite differences or finite elements (Particle-In-Cell, PIC). When simulating large-sized systems, this method encounters performance losses, both locally (intra-CPU) and globally (inter-CPU) because the amount of communicated data increases. Another approach, Particle-In-Fourier (PIF), takes advantage of the anisotropy of the field by representing it in Fourier space, where data is much smaller. One can scale up the system size keeping the number of particles per mode constant, and removing the subdomain decomposition, resulting in a better parallel scalability. The price to pay is the computation of transcendental functions at marker positions instead of low-order polynomials. We show that GPUs are well suited to such high arithmetic intensity. PIF versus PIC accuracy and performance comparisons are performed. Applications to physical cases involving mode coupling of high wave numbers are shown, for which the PIF approach performs faster than the PIC.