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} ' ' '0 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