% Listing 16.5 for Figure 16.5 of the book 
% "Principles of System Identification: Theory and Practice"

% Arun K. Tangirala
% December 06, 2015
% Freely distributed under no warranty

%% Frequency vector
Nfreq = 2000; winlen = 64;
wvec = (0:1/Nfreq:1-1/Nfreq)*2*pi - pi;

%% Windows and magnitudes of FRF
% Rectangular, Hann, Bartlett, Blackman, Flat top
w_rect = window(@rectwin,winlen); cg_rect = sum(w_rect);
w_hann = window(@hann,winlen); cg_hann = sum(w_hann);
w_bart = window(@bartlett,winlen); cg_bart = sum(w_bart);
w_black = window(@blackman,winlen); cg_black = sum(w_black);
w_flatt = window(@flattopwin,winlen); cg_flatt = sum(w_flatt);

figure; 
plot((0:winlen-1)/winlen,[w_hann w_bart w_black w_flatt]) 
hold on; plot((0:winlen-1)/winlen,w_rect,'k--')

%% Frequency-domain characteristics
Rectw_f = fft(w_rect/cg_rect,Nfreq);
Rectw_mag = abs(fftshift(Rectw_f));
dB_rect = 20*log10(Rectw_mag);
Hannw_f = fft(w_hann/cg_hann,Nfreq);
Hannw_mag = abs(fftshift(Hannw_f));
dB_Hann = 20*log10(Hannw_mag);
Bartw_f = fft(w_bart/cg_bart,Nfreq);
Bartw_mag = abs(fftshift(Bartw_f));
dB_Bart = 20*log10(Bartw_mag);
Blackw_f = fft(w_black/cg_black,Nfreq);
Blackw_mag = abs(fftshift(Blackw_f));
dB_Black = 20*log10(Blackw_mag);
Flattw_f = fft(w_flatt/cg_flatt,Nfreq);
Flattw_mag = abs(fftshift(Flattw_f));
dB_Flatt = 20*log10(Flattw_mag);

figure
plot(wvec/(2*pi/winlen),[dB_Hann dB_Bart dB_Black dB_Flatt])
hold on; plot(wvec/(2*pi/winlen),dB_rect,'k--');