%cleanSignal: Filter the raw received Rx signal.
% Usage: 
%       filtered_Rx_signal = filterSignal(Tx_raw_signal, raw_Rx_signal, f_Tx, f_sampling)

% See the following example for loading a database file and filter a
% specific sensor signal.

% START EXAMPLE
%% Select database file to read
file_name = '.\\data\\SARISTU\\50k\\A\\T0\\00\\calibA50k_00_11_pc.txt';
file_data = importfile(file_name);
file_data(1,:) = []; % removing Tx, Rxi, Rxj sensors numbers in file's first row
time = file_data(:,1);
file_data(:,1) = []; % removing time column
txr = file_data(:,1); % storing Tx signal
rxr = file_data(:,3); % Rx1 raw signal is in the 3rd column, Rx2 is in the 4th column.
                      % Check signal file to see Rx1 and Rx2 number.
clearvars file_data

%% Waveforms definition
fs = 2500000; % Sampling frequency, Hz, from file 10th row (SPS)
fc = 50e3; % Tx frequency signal, Hz

%% Calling filtering function
% rxr: raw signal to be filtered
% rxe: rxr filtered
rxe = filterSignal(txr, rxr, fc, fs);

%% Plotting
figure;
plot(time, rxr, time, rxe)
xlabel('time (s)')
ylabel('Voltage (V)')
legend('Raw signal', 'Filtered Signal')

% END EXAMPLE

%% --- Functions to perform the cleaning/filtering
function rxe = filterSignal(txr, rxr, fc, fs)
%% Options
opt.FrqScale=1; % Frequency scaling
opt.DownSmpl=1; % Downsampling factor
opt.TimeGate=1; % Enable/disable time gating
opt.tgWndExt=6; % Time window width multiplier
opt.tgWndTh=0.99; % Time window threshold
opt.HpFreq=fc/3;   % High pass filter cut-off frequency
opt.LpFreq=fc*3;   % Low pass filter cut-off frequency

[txe,fs] = preProcess(txr,fs,'Options',opt);
    
[val,idxP]=max(txe);
idxP = findPeakRef(txe,idxP,10);
    
wndLen = opt.tgWndExt*sum(txe>opt.tgWndTh*val);
gtd.wf = chebwin(wndLen,100)';
    
idxWnd = findPeakRef(gtd.wf,fix(wndLen/2),3);
gtd.tb = ((1:wndLen) - idxWnd)/fs;
gtd.tg = 0.*numel(rxr);
 
rxe = preProcess(rxr,fs,'Options',opt,'GateData',gtd);
       
end
   
   
%% --- fftf_LP.m
function y = fftf_LP(x,fs,lp)
N = numel(x);
k = ceil(N*lp/fs);

X = fft(x);
X(k+1:end-k+1) = 0;
y = ifft(X);
end

%% --- fftf_HP.m
function y = fftf_HP(x,fs,hp)
N = numel(x);
k = ceil(N*hp/fs);

X = fft(x);
X([1:k end-k+2:end]) = 0;
y = ifft(X,'symmetric');
end

%% --- findPeakRef.m
function idX = findPeakRef(mX,o,S)

    g = o+(-S:S);
    idX = g(find(mX(g)==max(mX(g)),1));
    
    while (idX>2) && (mX(idX-1) > mX(idX)), idX = idX-1; end;
    while (idX<numel(mX)-1) && (mX(idX+1) > mX(idX)), idX = idX+1; end;
    
    idX = idX + 0.5*(mX(idX-1)-mX(idX+1))/(mX(idX-1)-2*mX(idX)+mX(idX+1));
    
end


%% --- preProcess.m
function [sig,fs]=preProcess(sig,fs,varargin)

p = inputParser;
addOptional(p,'Options',setOptions); % Default options
addOptional(p,'Baseline',[]);        % Signal baseline
addOptional(p,'WarpMap',[]);         % Warping map
addOptional(p,'GateData',[]);        % Time gating data
parse(p,varargin{:});

% Store local options
opt = p.Results.Options;
bsl = p.Results.Baseline;
map = p.Results.WarpMap;
gtd = p.Results.GateData;
clear p;

% Force row vector format
if ~isrow(sig), sig = sig.'; end

% Bandpass filtering
sig = fftf_LP(sig, fs, opt.LpFreq);
sig = fftf_HP(sig, fs, opt.HpFreq);

% Time gating
if (opt.TimeGate) && ~isempty(gtd)
    t = (0:numel(sig)-1)/fs;
    msk = ones(1,numel(sig));
    for i=1:numel(gtd.tg)
        if gtd.tg(i)>numel(sig)-10, continue, end;
        tRef = gtd.tg(i)/fs;      
        msk = msk .* (1-interp1(gtd.tb,gtd.wf,t-tRef,'pchip',0));
    end
    
    sig = sig .* msk;
end

end

function data_file = importfile(filename, startRow, endRow)
%IMPORTFILE Import numeric data from a text file as a matrix.
%   CALIBA50K0001PC = IMPORTFILE(FILENAME) Reads data from text file
%   FILENAME for the default selection.
%
%   CALIBA50K0001PC = IMPORTFILE(FILENAME, STARTROW, ENDROW) Reads data
%   from rows STARTROW through ENDROW of text file FILENAME.
%
% Example:
%   calibA50k0001pc = importfile('calibA50k_00_01_pc.txt', [6,21], [6,2520]);
%
%    See also TEXTSCAN.

% Auto-generated by MATLAB on 2019/11/27 10:31:41

%% Initialize variables.
delimiter = '\t';
if nargin<=2
    startRow = [6,21];
    endRow = [6,2520];
end

%% Format for each line of text:
%   column1: double (%f)
%	column2: double (%f)
%   column3: double (%f)
%	column4: double (%f)
% For more information, see the TEXTSCAN documentation.
formatSpec = '%f%f%f%f%*s%*s%[^\n\r]';

%% Open the text file.
fileID = fopen(filename,'r');

%% Read columns of data according to the format.
% This call is based on the structure of the file used to generate this
% code. If an error occurs for a different file, try regenerating the code
% from the Import Tool.
dataArray = textscan(fileID, formatSpec, endRow(1)-startRow(1)+1, 'Delimiter', delimiter, 'TextType', 'string', 'EmptyValue', NaN, 'HeaderLines', startRow(1)-1, 'ReturnOnError', false, 'EndOfLine', '\r\n');
for block=2:length(startRow)
    frewind(fileID);
    dataArrayBlock = textscan(fileID, formatSpec, endRow(block)-startRow(block)+1, 'Delimiter', delimiter, 'TextType', 'string', 'EmptyValue', NaN, 'HeaderLines', startRow(block)-1, 'ReturnOnError', false, 'EndOfLine', '\r\n');
    for col=1:length(dataArray)
        dataArray{col} = [dataArray{col};dataArrayBlock{col}];
    end
end

%% Close the text file.
fclose(fileID);

%% Post processing for unimportable data.
% No unimportable data rules were applied during the import, so no post
% processing code is included. To generate code which works for
% unimportable data, select unimportable cells in a file and regenerate the
% script.

%% Create output variable
data_file = [dataArray{1:end-1}];

end

function out = addWave(txs,rxs,txw,rxw,rxb,fs,fc)
% addWave: Construction of the dataset to be used in locateDefect## functions.
% IN:   txs Tx sensor number
%       rxs Rx sensor number
%       txw signal trasmitted by txs, mx1
%       rxw signal received by sensor rxs, mx1
%       rxb baseline
%       fs sampling frequency, Hz
%       fc sampling signal, HZ. For numerical simulation fc=1/dt
%
% OUT:  out the data structure to pass to locateDefect## functions
%
% USAGE: See packData.m for a usage example.

out.txs = txs;
out.txw = txw;
out.rxs = rxs;
out.rxw = rxw;
out.rxb = rxb;
out.fs  = fs;
out.fc  = fc;
end

function opt = setOptions(varargin)
%setOptions define some options to pass to defectLocate## functions. See
% setOptions.m file for the list of possible settable options.
% 
% IN:   structure with new options. If not passed, default values are used
%
% OUT:  opt, the options data structure

p = inputParser;
addOptional(p,'FrqScale',1); % Frequency scaling
addOptional(p,'TimeGate',1); % Enable/disable time gating
addOptional(p,'tgWndExt',7); % Time window width multiplier
addOptional(p,'tgWndTh',.7); % Time window threshold
addOptional(p,'tgAdapt',1);  % Enable/disable adaptive time gating
addOptional(p,'RayTrcLv',1); % Set Raytracing depth
addOptional(p,'HpFreq',0);   % High pass filter cut-off frequency
addOptional(p,'LpFreq',0);   % Low pass filter cut-off frequency
parse(p,varargin{:});

opt = p.Results;
end








