function [varargout]=mat2wfdb(varargin)
%
% [xbit]=mat2wfdb(X,fname,Fs,bit_res,adu,info,gain,sg_name,baseline,isquant, isdigital)
%
% Convert data from a matlab array into Physionet WFDB format file.
%
% Input Paramater are:
%
% X       -(required)  NxM matrix of M signals with N samples each. The
%                      signals can be of type double.The signals are assumed to be
%                      in physical units already and will be converted to
%                      ADU.
% fname   -(required)  String where the the header (*.hea) and data (*.dat)
%          files will be saved (one single name for both, with no sufix).
% Fs      -(Optional)  1x1 sampling frequency in Hz (all signals must have
%          been sampled at the same frquency). Default is 1 Hz.
% bit_res -(Optional)  1xM (or Mx1):scalar determining the bit depth of the conversion for
%                      each signal.
%                      1x1 : If all the signals should have the same bit depth
%          Options are: 8,  16, and 32 ( all are signed types). 16 is the default.
% adu     -(Optional) Describes the physical units (default is 'mV').
%          Three input formats:
%            - String delimited by forward slashes (e.g. 'V/mV/mmHg'), with
%            M-1 slash characters
%            - Single string (e.g. 'V'), in which case all signals will 
%            have the same physical units.
%            - Cell array of strings, where the total units entered has to equal M 
%            (number of channels).
% info    -(Optional)  String that will be added to the comment section of the header file.
%           For multi-lined comments, use a cell array of strings. Each
%           cell will be output on a new line. Note that comments in the
%           header file are automatically prefixed with a pound symbol (#)
% gain    -(Required for digital only) Scalar or Mx1 array of floats indicating the difference in sample values 
%           that would be observed if a step of one physical unit occurred in the original 
%           analog signal. If the 'isdigital' field is 1, this field is mandatory. Otherwise,
%           this field is ignored if present. 
% baseline -(Required for digital only) Mx1 array of integers that specifies the sample value for each channel
%           corresponding to 0 physical units. Not to be confused with 'ADC zero' which 
%           is currently always taken and written as 0 in this function. If
%           the 'isdigital' field is 1, this field is mandatory. Otherwise,
%           this field is ignored if present. 
% sg_name -(Optional) Cell array of strings describing signal names.
%
% isquant   -(Optional) Logical value (default=0). Use this option if the
%           input signal is already quantitized and you want to remove round-off
%           error by mapping the original values to integers prior to fixed
%           point conversion. This field is only used for input physical
%           signals. If 'isdigital' is set to 1, this field is ignored.
%
% isdigital -(Optional) Logical value (default=0). Specifies whether the input signal is 
%            digital or physical (default). If it is digital, the signal values will be 
%            directly written to the file without scaling. If the signal is physical, 
%            the optimal gain and baseline will be calculated and used to digitize the signal
%            to write the WFDB file. This flag also decides the allowed
%            input combinations of the 'gain' and 'baseline' fields.
%            Digital signals must have both, and physical signals must have
%            neither (as the ideal values will be automatically calculated). 
%
% Ouput Parameter:
%
% xbit    -(Optional)  NxM the quantitized signals that written to file (possible
%          rescaled if no gain was provided at input). Useful for comparing
%          and estimating quatitization error with the input double signal X
%          (see examples below).
%
%
%  NOTE: The signals can have different amplitudes, they will all be scaled to
%  a reference gain, with the scaling factor saved in the *.hea file.
%
%Written by Ikaro Silva 2010
%Modified by Louis Mayaud 2011, Alistair Johson 2016
% Version 1.0
%
% Since 0.0.1
% See also wrsamp, wfdbdesc
%
%%%%%%%%%%  Example 1 %%%%%%%%%%%%
%
% display('***This example will write a  Ex1.dat and Ex1.hea file to your current directory!')
% s=input('Hit "ctrl + c" to quit or "Enter" to continue!');
%
% %Generate 3 different signals and convert them to signed 16 bit in WFDB format
% clear all;clc;close all
% N=1024;
% Fs=48000;
% tm=[0:1/Fs:(N-1)/Fs]';
% adu='V/mV/V';
% info='Example 1';
%
%
% %First signal a ramp with 2^16 unique levels and is set to (+-) 2^15 (Volts)
% %Thus the header file should have one quant step equal to (2^15-(-2^15))/(2^16) V.
% sig1=double(int16(linspace(-2^15,2^15,N)'));
%
% %Second signal is a sine wave with 2^8 unique levels and set to (+-) 1 (mV)
% %Thus the header file should one quant step equal a (1--1)/(2^8)  adu step
% sig2=double(int8(sin(2*pi*tm*1000).*(2^7)))./(2^7);
%
% %Third signal is a random binary signal set to to (+-) 1 (V) with DC (to be discarded)
% %Thus the header file should have one quant step equal a 1/(2^15) adu step.
% sig3=(rand(N,1) > 0.97)*2 -1 + 2^16;
%
% %Concatenate all signals and convert to WFDB format with default 16 bits (empty brackets)
% sig=[sig1 sig2 sig3];
% mat2wfdb(sig,'Ex1',Fs,[],adu,info)
%
% % %NOTE: If you have WFDB installed you can check the conversion by
% % %uncomenting and this section and running (notice that all signals are scaled
% % %to unit amplitude during conversion, with the header files keeping the gain info):
%
% % !rdsamp -r Ex1 > foo
% % x=dlmread('foo');
% % subplot(211)
% % plot(sig)
% % subplot(212)
% % plot(x(:,1),x(:,2));hold on;plot(x(:,1),x(:,3),'k');plot(x(:,1),x(:,4),'r')
%
%%%%%%%% End of Example 1%%%%%%%%%

%endOfHelp
machine_format='l'; % all wfdb formats are little endian except fmt 61 which this function does not support. Do NOT change this.
skip=0;

% Set default parameters
params={'x','fname','Fs','bit_res','adu','info','gain','sg_name','baseline','isquant', 'isdigital'};
Fs=1;
adu=[];
info=[];
isquant=0;
isdigital=0;
%Use cell array for baseline and gain in case of empty conditions
baseline=[];
gain=[];
sg_name=[];
x=[];
fname=[];
%Used to convert signal from double to appropiate type
bit_res = 16 ;
bit_res_suport=[8 16 32];

for i=1:nargin
    if(~isempty(varargin{i}))
        eval([params{i} '= varargin{i};']);
    end
end

disp(isdigital);
% Check valid gain and baseline combinations depending on whether the input is digital or physical.
if isdigital % digital input signal
    if (isempty(gain) || isempty(baseline))
        error('Input digital signals are directly written to files without scaling. Must also input gain and baseline for correct interpretation of written file.');   
    end
    if (~isempty(find(baseline>2147483647))||~isempty(find(baseline<-2147483648))) % baseline stored as int in wfdb library. 
        error('Baseline field must lie between 2^-31 and 2^31-1 for this WFDB version'); % Prevent bit overflow
    end
else % physical input signal
    if ( ~isempty(gain) || ~isempty(baseline)) % User inputs gain or baseline to map the physical to digital values.
        % Sorry, we cannot trust that they did it correctly... 
        warning('Input gain and baseline fields ignored for physical input signal. This function automatically calculates and applies the ideal values');
    end
end
    
switch bit_res % Write formats. 
    case 8
        fmt='80';
    case 16
        fmt='16';
    case 32
        fmt='32';
end

[N,M]=size(x);

if isempty(adu) % default unit: 'mV'
    adu=repmat({'mV'},[M 1]);
elseif iscell(adu) 
    % adu directly input as a cell array of strings
elseif ischar(adu)
    if ~isempty(strfind(adu,'/'))
        adu=regexp(adu,'/','split');
    else
        adu = repmat({adu},[M,1]);
    end
end

% ensure we have the right number of units
if numel(adu) ~= M
    error('adu:wrongNumberOfElements','adu cell array has incorrect number of elements');
end

if(isempty(gain))
    gain=cell(M,1); %Generate empty cells as default
elseif(length(gain)==1)
    gain=repmat(gain,[M 1]);
end
% ensure gain is a cell array
if isnumeric(gain)
    gain=num2cell(gain);
end

if(isempty(sg_name))
    sg_name=repmat({''},[M 1]);
end
if ~isempty(setdiff(bit_res,bit_res_suport))
    error(['Bit res should be one of: ' num2str(bit_res_suport)]);
end
if(isempty(baseline))
    baseline=cell(M,1); %Generate empty cells as default
elseif(length(baseline)==1)
    baseline=repmat(baseline,[M 1]);
end
% ensure baseline is a cell array
if isnumeric(baseline)
    baseline=num2cell(baseline);
end

if isempty(isquant)
    isquant = zeros(M,1);
elseif numel(isquant)==1
    isquant = repmat(isquant,[M,1]);
elseif numel(isquant)~=M
    error('isquant:wrongNumberOfElements','isquant  array has incorrect number of elements');
end


%Head record specification line
head_str=cell(M+1,1);
head_str(1)={[fname ' ' num2str(M) ' ' num2str(Fs) ' ' num2str(N)]};

switch bit_res % Allocate space for digital signals
    case 8
        y=uint8(zeros(N,M));
    case 16
        y=int16(zeros(N,M));
    case 32
        y=int32(zeros(N,M));
end

%Loop through all signals, digitizing them and generating lines in header file
for m=1:M
    nameArray = regexp(fname,'/','split');
    if ~isempty(nameArray)
        fname = nameArray{end};
    end
    
    [tmp_bit1,bit_gain,baseline_tmp,ck_sum]=quant(x(:,m), ...
        bit_res, gain{m}, baseline{m}, isquant(m), isdigital);
    
    y(:,m)=tmp_bit1;
    
    % Header file signal specification lines
    % Should we specify precision of num2str(gain)?
    head_str(m+1)={[fname '.dat ' fmt ' ' num2str(bit_gain) '(' ...
        num2str(baseline_tmp) ')/' adu{m} ' ' num2str(bit_res) ' 0 ' num2str(tmp_bit1(1)) ' ' num2str(ck_sum) ' 0 ' sg_name{m}]};
end
if(length(y)<1)
    error(['Converted data is empty. Exiting without saving file...']);
end

%Write *.dat file
fid = fopen([fname '.dat'],'wb',machine_format);
if(~fid)
    error(['Could not create data file for writing: ' fname]);
end

if (bit_res==8)
    count=fwrite(fid, y','uint8',skip,machine_format);
else
    count=fwrite(fid, y',['int' num2str(bit_res)],skip,machine_format);
end

if(~count)
    fclose(fid);
    error(['Could not data write to file: ' fname]);
end

fprintf(['Generated *.dat file: ' fname '\n']);
fclose(fid);

%Write *.hea file
fid = fopen([fname '.hea'],'w');
for m=1:M+1
    if(~fid)
        error(['Could not create header file for writing: ' fname]);
    end
    fprintf(fid,'%s\n',head_str{m});
end

if(~isempty(info))
    if ischar(info)
        fprintf(fid,'#%s\n',info);
    elseif iscell(info)
        for m=1:numel(info)
            fprintf(fid,'#%s\n',info{m});
        end
    end
end

if(nargout==1)
    varargout(1)={y};
end
fprintf(['Generated *.hea file: ' fname '\n']);
fclose(fid);

end

%%%End of Main %%%%


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%Helper function
function [y,adc_gain,baseline,check_sum]=quant(x, bit_res, gain, baseline, isquant, isdigital)

min_x=min(x(~isnan(x)));
max_x=max(x(~isnan(x)));
nan_ind=isnan(x);
rg=max_x-min_x;

if(isdigital) 
    % Digital input signal. Do not scale or shift the signal. The gain/baseline will only 
    % be used to write the header file to interpret the output wfdb record.
    if ((min_x < -2^(bit_res-1)+1) || (max_x > (2^(bit_res-1)-1 )))
        error(['Digital input signal exceeds allowed range of specified output format: {' num2str(-2^(bit_res-1)+1) ' < x < ' num2str(2^(bit_res-1)-1) '}']);
    end
    adc_gain=gain;
    y=x;
    
else
    % Physical input signal - calculate the gain and baseline to minimize
    % the detail loss during ADC conversion: y = gain*x + baseline. Ignore any input gain or baseline
    
    % Calculate the adc_gain, baseline, and map the signal to digital
    % Make sure baseline doesn't go beyond 4 byte integer range
    
    if rg==0 % Flatline signal. Manually set adc_gain or it will be infinite.
          baseline=-(2^(bit_res-1))+1; % Set baseline to minimum value of bit res
          if x(1)==0
            adc_gain=1; % Arbitrary gain=1 for all 0 input signal. All values stored as baseline. 
          else
            adc_gain=-baseline/x(1); % Set gain as inverse to store all values as exactly 0.  
          end
          
    else % Non flatline signal: adc_gain = (range of encoding / range of Data) -- remember 1 quant level is for storing NaN
        % Constraint - baseline must be stored as a 4 byte integer for the WFDB library. 
        if ((min_x>0) && (bit_res==32)) % All values are +ve, map with baseline = -2^31+1
            adc_gain=((2^32)-2)/max_x; 
            baseline=-2147483647;
            if isquant==0 % Only display warning message if not recalculating later
                disp('Due to baseline constraints, output precision may be slightly less than 32 bits for the positive channel.');
            end
        elseif((max_x<0) && (bit_res==32)) % All values are -ve, map with baseline = 2^31-1
            adc_gain=((2^32)-2)/abs(min_x);
            baseline=2147483647;
            if isquant==0 % Only display warning message if not recalculating later
                disp('Due to baseline constraints, output precision may be slightly less than 32 bits for the negative channel.');
            end
        else % Signal has both +ve and -ve values or fmt is not 32. Use full range of bits. 
            adc_gain=((2^bit_res)-2)/rg;
            baseline=round(-(2^(bit_res-1))+1-min_x*adc_gain);
            
        end
        if(isquant)
            % The (non flatline) input signal was already quantitized. Remove round-off error 
            % by setting the original values to integers prior to fixed point conversion
            
            xvalues=sort(unique(x(~isnan(x)))); % All the values of x
            incmin=min(diff(xvalues)); % An estimate of the smallest possible increment in the input signal
            quantlevels=rg/incmin; % The estimated number of quantization levels in the input signal
            
            % We want to map 1 increment to 1 digital unit. First make sure
            % the full increment range is less than the 2^N-2 increments able to be encoded
            % by the chosen bit resolution. The incmin estimate will always
            % be equal to or larger than the true incmin, so it won't
            % wrongly trigger errors in this validation step. 
            
            if (quantlevels>2^bit_res-2)
                if bit_res==32
                    disp(['The input signal has more quantization levels than 32 bits -1. ' ...
                        'Cannot directly map all input values to integers. Up to 1 bit of roundoff error may occur. Continuing...']);
                    calcquant=0; % Skip the integer matching and keep the old baseline/gain calculated. 
                else
                    error(['The input signal has more quantization levels than the chosen bit resolution. ' ...
                        'Please choose a higher resolution or remove the isquant option to allow up to 1 bit of roundoff error']);
                end
            else
                calcquant=1;
            end
            
            % Calculate gain+offset. Baseline must be stored as a 4 byte integer for the WFDB library.
            if calcquant
                if ((min_x>0) && (bit_res==32)) % 32 bit +ve quant mapping
                    adc_gain=1/incmin;
                    baseline=round(2147483647-adc_gain*max_x); % map max_x to 2^31-1. 
                    if (baseline<-2147483647) % Check if baseline goes below -2^31-1. If so, no quant. Recalculate gain and base. 
                        adc_gain=((2^32)-2)/max_x; 
                        baseline=-2147483647;
                        disp('Due to baseline constraints, the channel will not be quantized. Output precision may be less than 32 bits for the positive channel.');
                    end
                elseif((max_x<0) && (bit_res==32)) % 32 bit -ve quant mapping
                    adc_gain=1/incmin;
                    baseline=round(-2147483647-adc_gain*min_x); % map min_x to -2^31+1. 
                    if (baseline>2147483647) % Check if baseline goes above 2^31-1. If so, no quant. Recalculate gain and base. 
                        adc_gain=((2^32)-2)/abs(min_x); 
                        baseline=2147483647;
                        disp('Due to baseline constraints, channel will not be quantized. Output precision may be less than 32 bits for the negative channel.');
                    end
                else % Signal has both +ve and -ve values or fmt is not 32. Can use full range of bits.
                    adc_gain=1/incmin; % 1 digital unit corresponds to the smallest physical increment.
                    baseline=round(-(2^(bit_res-1))+1-min_x*adc_gain); % xmin still maps to ymin. xmax will not go beyond y limit, baseline should not go beyond 32 bit limits.  
                end
            end
        end
        % Check for 8 and 16 bit format 'baseline' field overflow. VERY
        % uncommon situation. Occurs if entire signal is +ve or -ve
        % with very high magnitude.
        if (baseline>2147483647) % Signal is all negative with large magnitude.
            warning('Large offset input channel entered. Output precision may be less than specified format for the negative channel.');
            baseline=2147483647; % Baseline is max int value, min_x maps to min bitres value. 
            adc_gain=(-(2^(bit_res-1))+1-baseline)/min_x;
        elseif (baseline<-2147483647) % Signal is all positive with large magnitude.
            warning('Large offset input channel entered. Output precision may be less than specified format for the positive channel.');
            baseline=-2147483647; % Baseline is min int value, max_x maps to max bitres value.
            adc_gain=(2^(bit_res-1)-1-baseline)/max_x;
        end
    end

    y=x*adc_gain+baseline;
    
end % signal is in digital range. adc_gain and baseline have been calculated. 

% Convert signals to appropriate integer type, and shift any WFDB NaN int values to 
% a higher value so that they will not be read as NaN's by WFDB
switch bit_res % WFDB will interpret the smallest value as nan. 
    case 8
        WFDBNAN=-128;
        y=int8(y); 
    case 16
        WFDBNAN=-32768;
        y=int16(y);
    case 32
        WFDBNAN=-2147483648;
        y=int32(y);
end
iswfdbnan=find(y==WFDBNAN); 
if(~isempty(iswfdbnan))
    y(iswfdbnan)=WFDBNAN+1;
end

%Set original NaNs to WFDBNAN
y(nan_ind)=WFDBNAN;

%Calculate the 16-bit signed checksum of all samples in the signal
check_sum=sum(y);
M=check_sum/(2^15);
if(M<0)
    check_sum=mod(check_sum,-2^15);
    if(~check_sum && abs(M)<1)
        check_sum=-2^15;
    elseif (mod(ceil(M),2))
        check_sum=2^15 + check_sum;
    end
else
    check_sum=mod(check_sum,2^15);
    if(mod(floor(M),2))
        check_sum=-2^15+check_sum;
    end
end

% Note that checksum must be calculated on actual digital samples for format 80,
% not the shifted ones. Therefore we only convert to real format now. 
if bit_res==8
    y=uint8(int16(y)+128); % Convert into unsigned for writing byte offset format. 
end

% Signal is ready to be written to dat file. 

end


function y=get_names(str,deli)

y={};
old=1;
ind=regexp(str,deli);
ind(end+1)=length(str)+1;
for i=1:length(ind)
    y(end+1)={str(old:ind(i)-1)};
    old=ind(i)+1;
end

end


