Matlab etude demonstrates spectrum amplitude error when frequency is not a product of DF by integer

unit_spectrum_et2.m

% Spectral analysis etudes
% This etude demonstrates spectrum amplitude error when frequency is not a
% product of df by integer
% - Try to change f1 variable  to see wrong amplitude on spectrum 
% - Try to change NFFT variable to see zero-padding effect.
% code by [email protected]

clear

% PARAMETERS
% -------------------------------------------------------------------------

% Sample rate (Hz)
Fs = 48e3;
% Time resolution (secs)
dt = 1/Fs;
% Number of signal samples
NSIG = 2^10;
% Number of FFT input samples (zero padded signal samples)
NFFT = 2^10;
% Actual frequency resolution (Hz)
df = Fs/NSIG;
% Interpolated frequency resolution (Hz)
dfi = Fs/NFFT;

% Central frequency (Hz)
f0 = 100*df;
% Frequency difference between neighbor tones (Hz)
f1 = 2.1*df;
% Number of tone signals
nf = 15;
% Tones frequencies (Hz)
fsig = (f0-(nf-1)/2*f1):f1:(f0+nf/2*f1);
% Tones amplitudes 0-Pk (Volts)
asig = ones(size(fsig)); % uniform
%asig = (1:numel(fsig)) / numel(fsig); % stairs

% GENERATOR
% -------------------------------------------------------------------------

% Generate tone signals
s = zeros(nf, NSIG);
t = (0:NSIG-1)*dt;
for i = 1:nf
    s(i,1:NSIG) = asig(i)*sin(2*pi*fsig(i)*t);
end

% PROCESSOR
% -------------------------------------------------------------------------

S = zeros(nf, NFFT/2-1);
for i = 1:nf
    % Calculate FFT magnitude
    x = zeros(1,NFFT);
    x(1:NSIG) = s(i,:);
    X = fft(x, NFFT);
    X = X(2:NFFT/2); % X(1) is DC
    X = abs(X);
    
    % Multiply by 2 for 0-Pk amplitude
    X = X/NSIG*2;
    
    % Store
    S(i,:) = X;
end

% RESULTS
% -------------------------------------------------------------------------

% Interpolated spectrum frequnecies (Hz)
f = (1:NFFT/2-1)*dfi;
f_khz = f*1e-3;
fsig_khz = fsig*1e-3;

% Find max for each spectrum
apk = zeros(1,nf);
fpk = zeros(1,nf);
for i = 1:nf
    [~, j] = max(S(i,:));
    apk(i) = S(i,j);
    fpk(i) = f_khz(j);
end
[fpk, sidx] = sort(fpk);
apk = apk(sidx);

figure;
plot(f_khz, S', '-');
hold on;
plot(f_khz, S', '.', 'Color', [0.5 0.5 0.5]);
plot(fpk, apk, 'LineWidth', 2, 'Color', [0 0 0]);
fb = max(fsig_khz) - min(fsig_khz);
fc = mean(fsig_khz);
xlim([fc-fb fc+fb]);
ylim([0 1.1*max(asig)]);
xlabel('kHz'); ylabel('V 0-Pk');
grid on; grid minor;