Matlab for Localized Peak Removal

        freq = 104.98;  % estimated peak frequency in Hz
        bw = 10;        % peak bandwidth estimate in Hz

        R = exp( - pi * bw / fs);            % pole radius
        z = R * exp(j * 2 * pi * freq / fs); % pole itself
        B = [1, -(z + conj(z)), z * conj(z)] % inverse filter numerator
        r = 0.9;     % zero/pole factor (notch isolation)
        A = B .* (r .^ [0 : length(B)-1]);   % inverse filter denominator

        residual = filter(B,A,bodyIR);       % apply inverse filter