% 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.
%
% CCK 2017-Feb-09 ported to Matlab (now runs in both Matlab & Octave).

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, postive, 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
stdout = 1; % Octave defines this constant by default.
fprintf(stdout, "%u chunks of %u floats, remainder %u loaded from file: %s \n", M, N, count, filename);

powerspec = coeffs2(1:N+1)/(2*M*N)^2;
   % DC thru Nyquist frequency.

fclose(datafile);