% WINDOW  Windowing tests by David Gibson
%
% This command plots a set of transforms of window functions
% and displays various parameters including bandwidth,
% noise-equivalent bandwidth, sidelobe height and coherent gain.
%
% The Matlab function FREQZ is used; but a second method, slower
% but more explicit, using FFT, is also documented. This method
% gives the same answers apart from a problem with a minus sign.
%
% The data for the Blackman-Harris windows does not agree with
% that given in Harris's paper. The BH window calculation uses
% COSINEWINDOW (based on BLACKMAN) and GENCOSWIN (based on GENCOSWIN)
% both written by David Gibson.
% 
% The latest release of MatLab includes a BLACKMANHARRIS function
% but my COSINEWINDOW is more generalised. I havent compared it with
% BLACKMANHARRIS because I do not have a valid upgrade license.


cla; clear ; format short compact
disp('-----------------------------------------------------------')

% ================== set parameters

n = 512;                    % length of window
k = 2^16;                   % number of sample points in FFT of window
sample_points_per_bin = k/n % print this just for info

% ================== calculate windows
    
w1 = boxcar(n);
w2 = blackman(n);         % [0.42    0.50    0.08    0]
w3 = cosinewindow(n, 'sym', [0.42323 0.49755 0.07922 0      ]);
w4 = cosinewindow(n, 'sym', [0.35875 0.48829 0.14128 0.01168]);
w5 = hamming(n);          % [0.54    0.46    0       0]

% OPTION:
% w5 = cosinewindow(n, 'sym', [0.35875 0.48829 0.14128 0.01174]);
% This w5 is a modified BH window, with a(4) adjusted so that 
% the coefficients sum to zero (a1 -a2 + a3 - a4) to give 18dB sidelobe rolloff

% OPTION:
% try a wideband window (peculiar result)...
% firstzero = 16; n2 = n / firstzero; % first zero for boxcar
% w5 = [ zeros((n-n2)/2,1); blackman(n2); zeros((n-n2)/2,1) ];

plots = 5;
ww = [w1 w2 w3 w4 w5];      % bundle the w's for convenience

[W1,f] = freqz(w1/sum(w1),1,k/2,2); W1 = W1 + 1E-10;  % to avoid log-of-zero errors
[W2,f] = freqz(w2/sum(w2),1,k/2,2);
[W3,f] = freqz(w3/sum(w3),1,k/2,2);
[W4,f] = freqz(w4/sum(w4),1,k/2,2);
[W5,f] = freqz(w5/sum(w5),1,k/2,2);
WW = [W1 W2 W3 W4 W5 ];

% ================== plot window transforms individually

filetype = '-dtiff';

% ================== window ===============
fig=1; figure(fig)
hh = plot( f*n/2,20*log10(abs(W1)));
grid;
set(gca,'LineWidth',1.5,'FontSize',12,'TickDir','out')
title('Rectangular Window Response','fontsize',12)
xlabel('Frequency (Bins)','fontsize',12)
ylabel('Normalized Magnitude (dB)','fontsize',12)

set (hh(1),'linewidth', 2)

ax = axis; 
axis([0 n/2 ax(3:4)])
axis([0 10 -90 0])

print(filetype, '', ['window-' num2str(fig)] )

% ================== window ===============
fig=2; figure(fig)
hh = plot( f*n/2,20*log10(abs(W3)));
grid;
set(gca,'LineWidth',1.5,'FontSize',12,'TickDir','out')
title('Blackman Harris 3-term Window Response','fontsize',12)
xlabel('Frequency (Bins)','fontsize',12)
ylabel('Normalized Magnitude (dB)','fontsize',12)

set (hh(1),'linewidth', 2)

ax = axis; 
axis([0 n/2 ax(3:4)])
axis([0 10 -90 0])

print(filetype, '', ['window-' num2str(fig)] )

break

% ================== plot window transforms all on one sheet

hh = plot( f*n/2,20*log10(abs([WW])));
grid;
legend('Rectangular','Blackman','Min 3-term B-H', 'Min 4-term B-H', ...
    'Hamming', 1) 
%   'Modified BH4', 1) 

title('Window Response') 
xlabel('Frequency (Bins)') 
ylabel('Normalized Magnitude (dB)') 

set (hh(1),'linewidth', 1)
set (hh(2),'linewidth', 1)
set (hh(3),'linewidth', 1.5)
set (hh(4),'linewidth', 2)
set (hh(5),'linewidth', 1)

ax = axis; 
axis([0 n/2 ax(3:4)])
axis([0 10 -120 0])

% ================== calculate bandwidth, sidelobe height, noise etc

% look for the 6dB bandwidth point.
% 6dB bandwidth is a more useful measure of resolving power than 3dB
% See p58 of Harris, 1978, 'On the Use of Windows for Harmonic Analysis 
% with the Discrete Fourier Transform', 
% Proc. IEEE, 66(1), Jan. 78, p51-83

[value index] = min(abs([WW]) > 1/sqrt(2));
bandwidth3dB = 2 * (index-1)*n/k
[value index] = min(abs([WW]) > 0.5);
bandwidth6dB = 2 * (index-1)*n/k

% look for position of first zero (or close to it)

[value index] = min(abs(WW) > 1E-5 );
%[value(5) index(5)] = min(abs(WW(:,5)) > 1E-3 );  % needed for hamming

zero = (index-1)*n/k

% Using the zero data, look for the maximum which occurs after that.
% This is the max sidelobe height.

for i=1:plots,
    [value2(i) index2(i)] = max(abs(WW(index(i):k/2, i)));
end
sidelobe = 20*log10(value2)  

% Now calculate the noise-equivalent bandwidth, using Harris'
% derivation, but noting that he omits a factor of N, possibly.

NEB = n * sum(ww.^2) ./ (sum(ww)).^2

% coherent gain, according to Harris, is

CG = sum(ww)/n

% scalloping loss is the magnitude of the point at a half-bin

SL = 20*log10( abs( WW(1 + sample_points_per_bin/2,:) ) )

% A suitable figure of merit is NEB / 3dB bandwidth

FOM = NEB ./ bandwidth3dB


% ================== results

% for the five plots above, the results are listed below
% with Harris's figures in []
%
%   Rect.     Blackman  BHarris3  BHarris4  Hamming
%
%bandwidth3dB =
%    0.8906    1.6563    1.6406    1.9063    1.3125
%   [0.89]    [1.68]    [1.66]    [1.90]    [1.30] 
%bandwidth6dB =
%    1.2188    2.3125    2.2813    2.6719    1.8281
%   [1.21]    [2.35]    [1.81]    [2.72]    [1.81] 
%zero =
%    1.0000    3.0000    3.0000    3.9844    2.0000
%sidelobe =
%  -13.2614  -58.1089  -70.2311  -92.0161  -42.6717
% [-13]     [-58]     [-67]     [-92]     [-43]
%NEB =
%    1.0000    1.7301    1.7118    2.0083    1.3647
%   [1]       [1.73]    [1.71]    [2.00]    [1.36]
%CG =
%    1.0000    0.4192    0.4224    0.3580    0.5391
%   [1.00]    [0.42]    [0.42]    [0.36]    [0.54]
%SL =
%    3.9224    1.0946    1.1244    0.8223    -1.7470
%   [3.92]    [1.1]     [1.13]    [0.83]     [1.78]
%FOM =
%    1.1228    1.0446    1.0434    1.0535    1.0398


% ================== repeat process, with different method
break


% instead of using FREQZ it is slightly more explicit to use FFT.
% The function we want is the ratio of the FFT of the window to the 
% FFT of an impulse.

% impulse = ifft(boxcar(n));
% That is just a convenience to generate [1 0 0 0 0 ...etc]
% in fact, lets not try to be so clever...

impulse = zeros(n,1); impulse(1) = 1;

% try inverse
%W1 = zeros(k,1); W1(1:100) = 1;
%w1 = ifft( fft((impulse),k) .* W1);

W1 = fft(w1/sum(w1),k) / fft((impulse),k) + 1E-10;
W2 = fft(w2/sum(w2),k) / fft((impulse),k) + 1E-10;
W3 = fft(w3/sum(w3),k) / fft((impulse),k) + 1E-10;
WW = [W1 W2 W3];
WW = WW(1:k/2);

% the problem is that this does not give the same answer. 
% WW (if I remember correctly) is the complex conjugate of what it should be.
% This does not affect the plots, but it is a mystery that ought to
% be solved.