Unverified Commit 9a0e85a9 authored by stephen scott tucker's avatar stephen scott tucker Committed by GitHub
Browse files

Generate filters with butter interface, raise warnings on NaN, inf

parent b9dc17dc
Loading
Loading
Loading
Loading
+30 −24
Original line number Diff line number Diff line
@@ -42,6 +42,7 @@ for ii=1:length(data)
        data2(ii) = AuxClass(data(ii));
    end
    y = data2(ii).GetDataTimeSeries();
    y2 = y;
    fs = data2(ii).GetTime();
    
    % convert t to fs
@@ -50,36 +51,41 @@ for ii=1:length(data)
        fs = 1/(fs(2)-fs(1));
    end
    
    % Check for NaN
    if max(max(isnan(y)))
       warning('Input to hmrR_BandpassFilt contains NaN values. Add hmrR_PreprocessIntensity_NAN to the processing stream.');
       return
    end
    % Check for finite values
    if max(max(isinf(y)))
       warning('Input to hmrR_BandpassFilt must be finite.');
       return
    end

    % Check that cutoff < nyquist
    if lpf / (fs / 2) > 1 || hpf / (fs / 2) > 1
        warning(['hmrR_BandpassFilt cutoff cannot exceed the folding frequency of the data with sample rate ', num2str(fs), ' hz.']);
        return
    end
    
    % low pass filter
    FilterType = 1;
    lpf_norm = lpf / (fs / 2);
    if lpf_norm > 0  % No lowpass if filter is 
        FilterOrder = 3;
    %[fa,fb]=butter(FilterOrder,lpf*2/fs);
    if FilterType==1 | FilterType==5
        [fb,fa] = MakeFilter(FilterType,FilterOrder,fs,lpf,'low');
    elseif FilterType==4
        %    [fb,fa] = MakeFilter(FilterType,FilterOrder,fs,lpf,'low',Filter_Rp,Filter_Rs);
    else
        %    [fb,fa] = MakeFilter(FilterType,FilterOrder,fs,lpf,'low',Filter_Rp);
        [z, p, k] = butter(FilterOrder, lpf_norm, 'low');
        [sos, g] = zp2sos(z, p, k);
        y2 = filtfilt(sos, g, double(y)); 
    end
    ylpf = filtfilt(fb,fa,double(y));
    
    % high pass filter
    FilterType = 1;
    hpf_norm = hpf / (fs / 2);
    if hpf_norm > 0
        FilterOrder = 5;
    if FilterType==1 | FilterType==5
        [fb,fa] = MakeFilter(FilterType,FilterOrder,fs,hpf,'high');
    elseif FilterType==4
        %    [fb,fa] = MakeFilter(FilterType,FilterOrder,fs,hpf,'high',Filter_Rp,Filter_Rs);
    else
        %    [fb,fa] = MakeFilter(FilterType,FilterOrder,fs,hpf,'high',Filter_Rp);
        [z, p, k] = butter(FilterOrder, hpf_norm, 'high');
        [sos, g] = zp2sos(z, p, k);
        y2 = filtfilt(sos, g, y2);
    end
    
    if FilterType~=5
        y2=filtfilt(fb,fa,ylpf);
    else
        y2 = ylpf;
    end
    data2(ii).SetDataTimeSeries(y2);
    
end