#-------------------------------------------------------------------------------
# @file      WindowFunction.m
# @brief     Octave/Matlab script for analysing different Windowing functions
# @author    Joern Henneberg
#            http://www.henneberg-online.org
#-------------------------------------------------------------------------------

clear all;         # clear all variables
hold off;   
figure(1); clf;    # clear figure
figure(2); clf;    # clear figure


#
# Calculate the Spectrum for the given window of size M
#
function spec = calc_spec(window, pad_factor, M, nfft)
    xw = [window, zeros(1, (pad_factor-1)*M)];      # create zero-padded window
    Xw = fft(xw);                                   # calc Fourier-Transform
    spec = 20*log10(abs(Xw));                       # Spectral magnitude in dB
    spec = spec - max(spec);                        # Normalize to 0 db max
    spec = max(spec,-350*ones(1,nfft));             # clip to -350 dB
endfunction


#
# Calculate a Rife-Vincent window of grade 4 and size N
#
function coeffs = rife_vincent4(N)
	# http://measure.feld.cvut.cz/usr/staff/haasz/publications/p32.pdf
	a0=35.0; a1=56.0; a2=28.0; a3=8.0; a4=1.0;    # Rife-Vincent4 coefficients
	coeffs = zeros(N, 1);
    for n = 1:N
        coeffs(n) = a0 - a1*cos(2*pi*n/N) + a2*cos(4*pi*n/N) - a3*cos(6*pi*n/N) + a4*cos(8*pi*n/N);
    endfor
	coeffs = coeffs / 128.0;                    # factor from Rife-Vincent formula
endfunction


# window size
M = 64;  

# calculate window coefficients
w_black    = blackman(M);          
w_blhar    = blackmanharris(M);    
w_hann     = hann(M);
w_hamm     = hamming(M);
w_kaiser3  = kaiser(M,3);
w_kaiser30 = kaiser(M,30);
w_rife4    = rife_vincent4(M);


zpf = 8;                            # zero-padding factor
nfft = zpf*M;                       # FFT size
fni = [0:1.0/nfft:1-1.0/nfft];      # Normalized frequency axis

# Calculate Spectrum for all windows
spec_black    = calc_spec(w_black'  , zpf, M, nfft);
spec_blhar    = calc_spec(w_blhar   , zpf, M, nfft);
spec_hann     = calc_spec(w_hann'   , zpf, M, nfft);
spec_hamm     = calc_spec(w_hamm'   , zpf, M, nfft);
spec_rife4    = calc_spec(w_rife4'  , zpf, M, nfft);
spec_kaiser3  = calc_spec(w_kaiser3', zpf, M, nfft);
spec_kaiser30 = calc_spec(w_kaiser30', zpf, M, nfft);


#
# Plot magnitude spectrum for well-performing windows
#
figure(1); 
plot(fni,spec_black,'b'); 
hold on;
plot(fni,spec_blhar,'r'); 
plot(fni,spec_rife4,'m'); 

axis([0,0.5,-250,10]); 
title('Comparison of well-perfoming Window Functions');
xlabel('Normalized Frequency'); 
ylabel('Magnitude (dB)');
legend('Blackman', 'Blackman-Harris', 'Rife-Vincent4');
grid on; 


#
# Plot magnitude spectrum for popular windows
#
figure(2); 
plot(fni,spec_hamm ,'m'); 
hold on;
plot(fni,spec_hann ,'g'); 
plot(fni,spec_black,'b'); 
plot(fni,spec_blhar,'r'); 

axis([0,0.5,-120,5]); 
title('Comparison of popular Window Functions');
xlabel('Normalized Frequency'); 
ylabel('Magnitude (dB)');
legend('Hamming', 'Hanning', 'Blackman', 'Blackman-Harris');
grid on; 


#
# end of file
#-------------------------------------------------------------------------------
