% FUNCTION ZEROCROSSING(FILENAME, [N_CHANNELS], [FRAME_LENGTH],
%                       [FRAME_SHIFT], [BLOCK_SIZE])
%
%   Compute zero-crossing rate frames of a multi-channel
%   waveform. Write out results in a file with same base as
%   FILENAME and extension '.zc%d' where %d is N_CHANNELS.
%
%   This function treats files in a block-per-block basis,
%   which allows for dealing with large files.
%
% INPUT:
%   filename:     string, name of the multi-channel waveform raw file.
%   n_channels:   integer, default 1. 
%   frame_length: integer, length of one frame in samples, default 512.
%   frame_shift:  integer, offset between two consecutive frames in
%                 samples, default floor(frame_length / 2).
%   block_size:   integer, number of frames to treat at a time,
%                 default 1000.
%
% FORMAT OF INPUT FILE:
%   Result of fwrite(fid, speech, 'int16', 'l') where speech is a
%   n_channels-by-n_samples matrix of signed 16-bit integers.
%   ('l' stands for little endian).
%
% FORMAT OF OUTPUT FILE:
%   Result of fwrite(fid, zc, 'float64', 'l') where zc is a 
%   n_channels-by-n_frames matrix of floats.
%   ('l' stands for little endian).
%
% ZERO-CROSSING RATE CALCULATION:
%   Per frame & per channel: number of zero-crossing is divided by
%   frame_length. Note that we weight the zero-crossing values (0
%   or 1) with hamming(frame_length - 1) before summing. This gives
%   more importance to zero-crossings in the middle of the frame.

function zerocrossing(filename, n_channels, frame_length, frame_shift, block_size)

% Default parameters
a = 2;    if nargin < a   n_channels = 1;   end
a = a+1;  if nargin < a   frame_length = 512; end
a = a+1;  if nargin < a   frame_shift = floor(frame_length / 2); end
a = a+1;  if nargin < a   block_size = 1000; end

% Open input file & get size
fid_in = fopen(filename, 'r');
if fid_in < 0
  error(sprintf ...
	('zerocrossing(): could not open waveform file %s.\n', filename));
end

fseek(fid_in, 0, 'eof');
n_samples = ftell(fid_in) / n_channels / 2;
fseek(fid_in, 0, 'bof');

if n_samples < frame_length
  error('zerocrossing(): not enough samples for even one frame!');
end

% Open output file
[p,n] = fileparts(filename);
filename_out = fullfile(p, [n '.zc' sprintf('%d', n_channels)]);
fid_out = fopen(filename_out, 'w');
if fid_out < 0
  error(sprintf ...
	('zerocrossing(): could not open output file %s.\n', filename_out));
end

% Proceed through data block per block (this allows for dealing
% with large files)

n_frames = 1 + floor((n_samples - frame_length) / frame_shift);
disp(sprintf('%d frames to treat', n_frames));

m = repmat(1 + (0:block_size-1) * frame_shift, frame_length , 1);
m = m + repmat((0:frame_length-1)', 1, block_size);

w = repmat(hamming(frame_length - 1), 1, block_size);

last_frame_shift = [];

for a = 1:block_size:n_frames
  
  n = min(n_frames, a+block_size-1) - a + 1;
  disp(sprintf('treating frames %d to %d', a, a+n-1));
  
  % Read block of data
  aux  = frame_length + (n - 1 - (1-isempty(last_frame_shift)) ) * frame_shift;
  [data, count] = fread(fid_in, [n_channels aux], 'int16', 'l');
  if count < aux * 8
    error('zerocrossing(): reached EOF before expected');
  end
  data = [last_frame_shift, data];
  
  % For each channel, compute zero-crossing rate of each frame
  for b = 1:n_channels
    x = data(b, :) / 32768; 
    y = sign(x(m(:, 1:n)));    % frame_length-by-n  
    y = .5 * abs(diff(y, 1, 1));
    y = y .* w(:, 1:n);
    zc(b,1:n) = sum(y, 1) / frame_length;
  end
  
  % Append results to the output file
  fwrite(fid_out, zc(:,1:n), 'float64', 'l');
  
  % Solve overlap between blocks
  last_frame_shift = data(:, end-frame_shift+1:end);
end

% Close files
fclose(fid_in);
fclose(fid_out);