Matlab etude demonstrates window function effect on multi-tone signal

unit_spectrum_et4.m

% Spectral analysis etudes
% This etude demonstrates window function effect on multi-tone signal
% - Try to change wtype variable to see the 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^14;
% 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 = 5*df;
% Number of tone components in signal
nf = 5;
% 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

% Window function
% 'rectwin' 'hamming' 'hann' 'flattopwin' 'blackmanharris'
wtype = 'hamming';

% 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
% -------------------------------------------------------------------------

% Calculate win. coef.
w = window(wtype, NSIG)';
w = w/mean(w);

% Calculate FFT magnitude
x = [sum(s,1).*w zeros(1,NFFT-NSIG)];
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;

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

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

% plot spectrum
figure; stem(f_khz, X, '.');
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;