High-frequency spur after performing an FFT in MATLAB

I have a modulated signal, and now I want to perform an FFT. However, I am getting a spur at a high frequency, which should not be there (and if it should, I have no clue as to why).

Lvl=[0.5,0.9,0.5,0.5,0.1,0.1,0.9,0.5];
fa=60;  %the frequency of the parasitic source in hertz
np=2;   %number of periods per bit

kl=length(Lvl);
t=0:0.01*np/fa:np*kl/fa;

Sig=sin(2*pi*fa*t);

for n=1:1:101
    Sig(n)=Sig(n)*Lvl(1);
end 

for n=102:1:201   
    Sig(n)=Sig(n)*Lvl(2);
end

for n=202:1:301
    Sig(n)=Sig(n)*Lvl(3);
end

for n=302:1:401
    Sig(n)=Sig(n)*Lvl(4);
end

for n=402:1:501
    Sig(n)=Sig(n)*Lvl(5);
end

for n=502:1:601
    Sig(n)=Sig(n)*Lvl(6);
end

for n=602:1:701
    Sig(n)=Sig(n)*Lvl(7);
end

for n=702:1:801
    Sig(n)=Sig(n)*Lvl(8);
end

plot(t,Sig)

%FFT
y = fft(Sig);     
f = (0:length(y)-1)*(1/(0.01*np/fa))/length(y);

plot(f,abs(y))
title('Magnitude')

I'm expecting just a spike at 60Hz with spurs around it, but instead I'm getting that and a large spike at almost 3kHz with spurs around it.

This peak at almost 3 kHz should be there, since the fft of a real is signal symmetric around the nyquist frequency (actually complex conjugate). The nyquist frequency is half the samping frequency, in your case sampling is done at 3000 Hz, thus the nyquist frequency is 1500 Hz. If you look closer at the peak, you will see that it is at 2940 Hz (which is 3000-60 Hz), due to the fact that the fft is mirrored around 1500 Hz. There are plenty of sources that explain why this is a property of the Fourier transform (e.g. here).

The actual Fourier transform would be mirrored around the zero frequency, but the fft gives you the fast Fourier transform, which is mirrored around the nyquist frequency. You can use fftshift to center the spectrum around the zero frequency.

I took the liberty to shorten your code, by avoiding repetition of several for-loops, and added the fftshift. Since your signal is real, you can also choose to show only one side of the fft, but I'll leave that up to you.

Lvl=[0.5,0.9,0.5,0.5,0.1,0.1,0.9,0.5];
fa=60;  % the frequency of the parasitic source in hertz
np=2;   % number of periods per bit

kl = length(Lvl);
dt = 0.01*np/fa; % time step 
Tend = np*kl/fa - dt; % time span 
t = 0:dt:Tend; % time vector
N = length(t); % number samples

Sig=sin(2*pi*fa*t);
for n = 1:kl
    ids = (1:100) + (n-1)*100; 
    Sig(ids) = Sig(ids).*Lvl(n);
end

% FFT
Y = fft(Sig);     
fv = (0:N-1)/(N*dt);     % frequency range
% FFT shift: 
Y_shift = fftshift(Y);
fv_shift = (-N/2:N/2-1)/(N*dt); % zero centered frequency vector

% Plot
figure(1); clf
subplot(311)
plot(t,Sig)
title('Signal')

subplot(312)
plot(fv,abs(Y))
title('FFT Magnitude')

subplot(313)
plot(fv_shift,abs(Y_shift))
title('FFT Magnitude zero shift')