% spectrum.m
function [powerspec,frequency] = spectrum(filename, N, dt)

% Power spectral estimation by partitioning and windowing. The program is 
% able to analyze data sets larger than memory by processing a binary stream.
% Big-Endian ieee floating point numbers (32 bits) are assumed.
% Binary file given by filename is read in chunks of N floats,
% assumed big-Endian. Power spectrum is built 2N samples at a time.
% The sampling period is optionally given by dt.

if (exist('dt') ~= 1) dt=1.0;
end

window = ( 1 - cos( 2*pi*(1:(2*N))/(2*N) ) ) / 2;  % Hann window

frequency = (0:N) / (2*N*dt);  % Convenient frequency axis
   % (cycles per unit time, with same time units as dt).
   % Note that only the DC, postivie, and Nyquist frequencies are represented here;
   % negative frequencies are redundant.

coeffs2 = zeros(1,2*N);  % Create empty array for squared Fourier coefficients

datafile = fopen(filename,"r");
% (We now need a do-until loop, which exists in Octave and has finally
% been added to MATLAB, but unfortunately with a different syntax. The
% following code, while slightly awkward, should work in both languages.)
M=0;
while (1)
   [chunk, count] = fread(datafile, [1,N], "float", "ieee-be");
   if (count ~= N) break; end;
   if (exist('lastchunk') == 1)
      coeffs2 = coeffs2 .+ abs( fft(window.*[lastchunk,chunk]) ).^2;   
         % matlab lacks +=.   
   end
   lastchunk = chunk;
   M=M+1;
end
fprintf(stdout, "%u chunks of %u floats, remainder %u loaded from file: %s \n", M, N, count, filename);

powerspec = [ coeffs2(1), 2*coeffs2(2:N), coeffs2(N+1) ]/(2*M*N)^2;
   % 1st term: DC
   % 2nd term: positive frequencies, doubling since an equal contribution
   %     is due from N+2:2N.
   % 3rd term: Nyquist frequency.
   % (cf. NR eqs. 13.4.5.)

fclose(datafile);