Best way to remove noise but preserve transients other than LPF

Question

Here's some code.

Fs = 10000;
t = linspace(0,10,10*Fs);
f = 100;
i  = zeros(size(t));
imp = [0 0.12 -0.2 0.1 1 -0.7 -0.2 0.1 0];
impTime = 0.2;
nSamp = impTime*Fs;
impFull = linspace(0,impTime,nSamp);
impUse = interp1(linspace(0,impTime,length(imp)),imp,impFull,'spline');
plot(impFull,impUse)
rate = 2;
nCycles = 10/rate;
cycleLen = rateFs;
cycle=zeros(1,cycleLen);
cycle(1:impTimeFs) = impUse;
cycles = repmat(cycle,1,10/rate);
x = 0.2sin(2pif/4t);
x1 = 0.1sin(2pift/3);
n = 0.3*randn(size(x));
y =x+x1+n;
y = y+2*cycles;
plot(t,y);

How would you go about removing the noise whilst preserving the transients, other than low pass filtering? And conversely, how would you remove the transients but preserve the noise?

Answer 1

In this case since the transients appear to be at a consistent rate (if that's the case) then a simple harmonic bandpass filter would be a good consideration. The filter would be a multi-harmonic resonator with the first harmonic at the periodic rate of the transients and repeat at multiple higher harmonics (the more that are included the more an impulse like response is preserved). If we have the opportunity for the sampling rate to be a multiple of the transient rate, then the filter is easily constructed as an IIR filter with the transfer function given as:

$$H(z) = \frac{1-\alpha}{1 - \alpha z^{-N}}$$

Where $N/2$ corresponds to the number of Harmonics to include between DC and half the sampling rate. $\alpha$ is related to the filter Q with a value between 0 and 1; the closer to 1, the higher the Q (the tighter the resonance), but also the more digital precision is required in implementation.

Below demonstrates the implementation with $N=10$ and $\alpha=0.999$.

What is compelling about this approach is its simplicity- it is a delay line with feedback:

For cases when the sampling rate is not an integer multiple of the transient rate, the same structure can be used but the integer sample delay $z^{-N}$ would be replaced with a fractional delay filter.

To remove the transients and preserve the noise, a harmonic nulling filter can be used which has the following transfer and frequency response, here demonstrated with $N=10$ and $\alpha=0.9$):

$$H(z) = \frac{1+\alpha}{2}\frac{1-z^{-N}}{1-\alpha z^{-N}}$$

And the block diagram showing the implementation of the harmonic nulling filter is shown below:

Answer 2

You have a few options to do what you want. If your primary goal is to denoise the signal but keep the transients then you can use wavelet denoising.

This is what the results look like:

If you want to remove the transients but keep the noise then you can simply subtract the denoised signal from the original signal to get a noisy signal with no transients:

But since this requires you to first denoise the signal and the subtract it from the original, you could directly remove the transients by high pass filtering:

Here is the code:

clear;
Fs = 10000;
t = linspace(0,10,10*Fs);
f = 100;
i  = zeros(size(t));
imp = [0 0.12 -0.2 0.1 1 -0.7 -0.2 0.1 0];
impTime = 0.2;
nSamp = impTime*Fs;
impFull = linspace(0,impTime,nSamp);
impUse = interp1(linspace(0,impTime,length(imp)),imp,impFull,'spline');
figure(1);
plot(impFull,impUse)
rate = 2;
nCycles = 10/rate;
cycleLen = rateFs;
cycle=zeros(1,cycleLen);
cycle(1:impTimeFs) = impUse;
cycles = repmat(cycle,1,10/rate);
x = 0.2sin(2pif/4t);
x1 = 0.1sin(2pift/3);
n = 0.3*randn(size(x));
y =x+x1+n;
y = y+2*cycles;
threshold = 2; % Adjust this threshold in INTEGER amounts only
wavelet = 'db4'; % You can choose a different Wavelet
denoised_y = wdenoise(y, threshold, 'Wavelet', wavelet);
figure(2);
subplot(2, 1, 1);
plot(t, y);
title('Original Noisy Signal');
subplot(2, 1, 2);
plot(t, denoised_y);
title('Denoised Signal');
% Transient Removal (Indirectly High-Pass Filtering)
transient_removed_y = y - denoised_y;
figure(3);
subplot(2, 1, 1);
plot(t, y);
title('Original Noisy Signal');
subplot(2, 1, 2);
plot(t, transient_removed_y);
title('Transient Removed Signal');
% Transient Removal (Directly High-Pass Filtering)
cutoff_frequency = 17/(Fs/2);
[b, a] = butter(6, cutoff_frequency, 'high'); % Not necessary to use Butterworth filter
butter_transient_removed_y = filtfilt(b, a, y);
figure(4);
subplot(2, 1, 1);
plot(t, y);
title('Original Noisy Signal');
subplot(2, 1, 2);
plot(t, butter_transient_removed_y);
title('Transient Removed Signal (Butterworth)');
MAE = compute_mae(butter_transient_removed_y, transient_removed_y, 100000); % Maximum Absolute Error between the two methods
fprintf("MAE between the transient removal methods is %.8f \n", MAE)
function mae = compute_mae (y1 , y2 ,N)
    for el =1:1: N
         mae = (abs (y1(el)-y2(el)));
    end
    mae = max(mae);
end

I have also added a measure for computing the difference between the two transient removal methods. You can play with the different parameters as you see fit to achieve even better results.

If you are looking for some other alternatives then this answer might also help!