To simulate seismic wavefields with a frequency-domain wave equation, conventional numerical methods must solve the equation sequentially to obtain the wavefields for different frequencies. The monofrequency equation has the form of a Helmholtz equation. When solving the Helmholtz equation for seismic wavefields with multiple frequencies, a physics-informed neural network (PINN) can be used. However, the PINN suffers from the problem of spectral bias when approximating high-frequency components. We propose to simulate seismic multi-frequency wavefields using a PINN with an embedded Fourier feature. The input to the Fourier feature PINN for simulating multi-frequency wavefields is four-dimensional, namely the horizontal and vertical spatial coordinates of the model, the horizontal position of the source, and the frequency, and the output is multi-frequency wavefields at arbitrary source positions. While an effective Fourier feature initialization strategy can lead to optimal convergence in training this network, the Fourier feature PINN simulates multi-frequency wavefields with reasonable efficiency and accuracy.