% [ 2006-CHN-USS/MAIN.m ]
%
% This is an online implementation of the RSE-USS approach,
% published in ASRU 2005 (Lathoud, Magimai-Doss, Bourlard).
%
% The main function is "uss_filter.m".
%
% It reads a waveform block by block (e.g. 1-second block),
% and for each block:
%
% - optionally, CHannel Normalization is performed (CHN).
%
% - optionally, a Rayleigh + Shifted Erlang model is fitted (RSE),
%   and from the noise level parameter in the RSE model,
%   Unsupervised Spectral Subtraction is applied (USS).
%
% - optionally, GMN-USS is applied. This is a simplified 
%   version of RSE-USS, where the noise level is estimated
%   as the geometric mean of all observed magnitude values.
%
% It can have a most 3 outputs:
% - HTK file containing magnitude spectrum.
% - HTK file containing filterbank values.
% - HTK file containing MFCC values.
%
% For filterbanks and MFCCs, you may want to select values
% that are appropriate to your data:
% - par.mfcc_band_hz
% - par.mfcc_B
% - par.mfcc_N
% - par.mfcc_L
%
% By Guillaume Lathoud 2006
% lathoud@idiap.ch


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% First example: Aurora file
% ( raw, big endian, 8000 Hz )

par = [];

%par.input_filename = 'DATA/Aurora-A-N1SNR15-FCJ_Z2319Z9A.08';
par.input_filename = 'DATA/Aurora-C-N1SNR15-FCJ_Z2319Z9A.08';
par.input_format   = 'raw-be';
par.fs             = 8000;


par.htkout_magspectrum_filename = '/tmp/out.magspectrum_htk';
% par.htkout_magspectrum_filename = [];  % No output

par.htkout_fbank_filename       = '/tmp/out.fbank_htk';
% par.htkout_fbank_filename       = [];  % No output

par.htkout_mfcc_filename        = '/tmp/out.mfcc_htk';
% par.htkout_mfcc_filename        = [];  % No output

par.frameshift_sec        = 0.010;
par.framelength_sec       = 0.032;

par.block_size_sec = 1.0;

par.do_chn = 1;

par.do_rse_uss = 1;
par.do_gmn_uss = 0;

%par.mfcc_band_hz = [ -Inf  +Inf ];  % fullband
%par.mfcc_band_hz = [ 300 3400 ];  % telephone band
par.mfcc_band_hz = [ 180 3800 ]; % custom

par.debug   = 0;
par.verbose = 1;

dbstop error

uss_filter( par );

% Show the result

mag = my_htk_read( '/tmp/out.magspectrum_htk' );
figure; imagesc( log( mag ) ); title( 'log mag' );

fbank = my_htk_read( '/tmp/out.fbank_htk' );
figure; imagesc( log( fbank ) ); title( 'log fbank' );

mfcc = my_htk_read( '/tmp/out.mfcc_htk' );
figure; imagesc( mfcc ); title( 'mfcc' );


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Second example: meeting room file
% ( wav, 16000 Hz )

par = [];

%par.input_filename = 'DATA/av163-meetingroom-distant.wav';
par.input_filename = 'DATA/clean-NU-1028.streetaddr.with-zeroes.wav';  % Sanity check
par.input_format   = 'wav';
% par.fs             % Not necessary: automatically read

par.htkout_magspectrum_filename = '/tmp/out.magspectrum_htk';
par.htkout_fbank_filename       = '/tmp/out.fbank_htk';
par.htkout_mfcc_filename        = '/tmp/out.mfcc_htk';
par.frameshift_sec        = 0.010;
par.framelength_sec       = 0.032;

par.block_size_sec = 1.0;

par.do_chn = 1;

par.do_rse_uss = 1;
par.do_gmn_uss = 0;

par.mfcc_band_hz = [ -Inf  +Inf ];  % fullband
%par.mfcc_band_hz = [ 300 3400 ];  % telephone band
%par.mfcc_band_hz = [ 180 3800 ]; % custom

par.debug   = 0;
par.verbose = 1;

dbstop error

uss_filter( par );

% Show the result

mag = my_htk_read( '/tmp/out.magspectrum_htk' );
figure; imagesc( log( mag ) ); title( 'log mag' );

fbank = my_htk_read( '/tmp/out.fbank_htk' );
figure; imagesc( log( fbank ) ); title( 'log fbank' );

mfcc = my_htk_read( '/tmp/out.mfcc_htk' );
figure; imagesc( mfcc ); title( 'mfcc' );


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Third example: meeting room file
% ( NIST SPHERE format, 16000 Hz )

par = [];

par.input_filename = 'DATA/av163-meetingroom-distant.sph';
%par.input_filename  = 'DATA/clean-NU-1028.streetaddr.sph';

par.input_format   = 'sphere';

% par.fs             % Not necessary: automatically read

par.htkout_magspectrum_filename = '/tmp/out.magspectrum_htk';
par.htkout_fbank_filename       = '/tmp/out.fbank_htk';
par.htkout_mfcc_filename        = '/tmp/out.mfcc_htk';
par.frameshift_sec        = 0.010;
par.framelength_sec       = 0.032;

par.block_size_sec = 1.0;

par.do_chn = 1;

par.do_rse_uss = 1;
par.do_gmn_uss = 0;

%par.mfcc_band_hz = [ -Inf  +Inf ];  % fullband
%par.mfcc_band_hz = [ 300 3400 ];  % telephone band
par.mfcc_band_hz = [ 180 +Inf ]; % custom

par.debug   = 0;
par.verbose = 1;

dbstop error

uss_filter( par );

% Show the result

mag = my_htk_read( '/tmp/out.magspectrum_htk' );
figure; imagesc( log( mag ) ); title( 'log mag' );

fbank = my_htk_read( '/tmp/out.fbank_htk' );
figure; imagesc( log( fbank ) ); title( 'log fbank' );

mfcc = my_htk_read( '/tmp/out.mfcc_htk' );
figure; imagesc( mfcc ); title( 'mfcc' );
