## Spectrum analysis of the input file
## John A. Phillips, 21 August 2005

# Some constants
# window = hanning(Lfft) is built-in below;
Koverlap = 0.5;
Nwindows = 4;
FFToption = 0;
videowin2 = 1;
fname = "test.wav"
fdir = "/home/john/projects/resample-test/";

# Read the signal from the file.  The signal is returned as a 
# real value in the range -1.0 to +1.0 (approximately) so 
# there's no need to normalize.
[signal, Fsample, Nbits] = wavread(strcat(fdir,fname));
Nsamples = rows(signal);
Nchannels = columns(signal);

# Fit the required number of overlapping windows into the number
# of samples we have (FFToption == 0). 
Lwindow = floor(Nsamples/(Nwindows-Koverlap*Nwindows+Koverlap));
Lfft = Lwindow;

# If we need to have the length of the FFT in the form 2^N we 
# have two options: 1. reduce Lwindow and Lfft to the next power 
# of 2 below their calculated sizes, possibly using more windows
# as a result; or 2. make Lfft the next power of 2 above Lwindow 
# and pad the FFT with zeroes (which automatically happens later). 
if FFToption == 1
  Lwindow = 2^(nextpow2(Lwindow+1)-1);
  Lfft = Lwindow;
endif
if FFToption == 2
  Lfft = 2^(nextpow2(Lwindow));
endif

# Compute the window offsets into the signal file
overlap = ceil(Koverlap*Lwindow);
offsets = 1 : Lwindow-overlap : Nsamples-Lwindow+1;

# What's the total source power of the bit we will use?
# s = offsets(1);
# e = offsets(length(offsets)) + Lwindow - 1;
# source_power_dB = 10.0*log10(max(realmin,sumsq(signal(s:e,:),1)/(e+1-s)))

# Create the power-normalized spectral window
w = hanning(Lwindow);
w = w / norm(w);
window = zeros(Lwindow,Nchannels);
for i = 1:Nchannels window(:,i) = w; endfor

# Average over all of the windows
power = zeros(Lfft,Nchannels);
for i = 1:length(offsets)
  wsignal = signal(offsets(i):offsets(i)+Lwindow-1,:); 
  wsignal = fft(postpad(wsignal.*window,Lfft,0),Lfft,1); 
  power = power + real(wsignal.*conj(wsignal));
endfor
power = power / length(offsets);

# Throw away the duplicate half of the FFT, compensate for having 
# thrown away half of the energy and compensate for the length of 
# the FFT
Npoints = ceil((Lfft+1)/2); 
y = 2.0/Lfft * power(1:Npoints,:); 

# However in the simple energy compensation above, the zero-
# frequency point was not duplicated in both halves of the FFT so
# we must un-compensate it ...
y(1,:) = y(1,:)/2.0; 

# ... and if Lfft is even the Nyquist point was also not duplicated 
# in both halves of the FFT.  If so we have to un-compensate it too
if (rem(Lfft,2) == 0)
  y(Npoints,:) = y(Npoints,:)/2.0; 
end 

# Check that the total FFT power still equals the signal power, 
# but note that transient signals in the file can spoil this check
# FFT_power_dB = 10.0 * log10(max(realmin,sum(y,1)))

# Increase the video bandwidth if the half-video-bandwidth window 
# videowin2 is set to > 0. This compensates for the spreading to 
# make spectral peak levels correct although noise levels increase.
if (videowin2 > 0)
  video = zeros(Npoints+2*videowin2,Nchannels);
  for i = 1:2*videowin2+1
    video(i:i+Npoints-1,:) = video(i:i+Npoints-1,:) + y;
  end
  y = video(1+videowin2:Npoints+videowin2,:);
endif

# Create the frequency vector.  Points in the original FFT start 
# at zero frequency and are spaced by Fsample/Lfft.  Now we have 
# only Npoints at the same start frequency and spacing
f = (0:Npoints-1) * Fsample/Lfft; 

# Generate the plot, title and labels
title(""); 
xlabel("Frequency (Hz)"); 
ylabel("Power (dBFS)"); 
grid("on");
for i = 1:Nchannels
  if (Nchannels == 1) 
    format = strcat(";",fname,";"); 
  elseif (Nchannels == 2)
    if (i == 1) format = strcat(";",fname," left ch.;"); endif
    if (i == 2) format = strcat(";",fname," right ch.;"); endif
  else
    format = strcat(";",fname," ch.",int2str(i),";");
  endif
  plot(f, 10.0 * log10(max(realmin,y(:,i))),format); 
  hold on
endfor
hold off