Which filtering approach to use?

I am dealing with a transient signal like this:

It is essentially zero in the beginning and it is important to keep this baseline.

I would like to get rid of the low frequency part in it, say to apply a 10kHz high-pass filter. There are commercial codes, which can do such filtering, but I am struggling to find a proper Matlab solution.

The best I found so far is a sos-based filtfilt approach, but this still results in something like this:

Thus it does not preserve the baseline, but adds artifical oscillations.

Any help or ideas for alternatives are welcome!

Kind regards, Markus

[Chad Greene]

What do you mean by preserving the baseline? It looks like the second signal does pass through (nearly) y=0 at x=0. If you want to preserve the overall shape, maybe you're looking for a bandstop filter?

[Markus Sause]

The "baseline" from 0 ... 30µs is approximately zero. Applying a filter to this section should return the same zero. So the oscillation in the filtered signal is purely artificial, likely because of the mismatch of periodic nature of the signal.

Answer 1

[Star Strider]

You cannot filter out the low-frequency part of your signal with a highpass filter and retain the constant, 0 Hz, section for the first 40 µs. That is a logical (and mathematical) impossibility. You must use a lowpass filter if you want to retain the constant section.

What do you want to do with your signal? What information in it do you want to keep, and what do you want to eliminate?

Attaching a .mat file with your time and amplitude vectors would help.

[Markus Sause]

I do agree partially. It is theoretically possible and is already technically implemented. But I have never seen another software besides this commerical one being able to do that, so I was wondering, which algorithm would be capable of doing that. For comparison, here are screenshots of what I would like to achieve (signal was normalized as this software cannot deal with the floating point precision)

Before filtering:

After filtering:

According to the software settings, this is done by a "Butterworth, 6th order, 100kHz". I hardly think so and assume there is more than that. If someone wants to give it a try, attached is the data file.

Sorry for the aliasing!

[Star Strider]

I downloaded the ‘simwaves.txt’ file and will explore it. It could take time to duplicate the software result.

If you could also upload the ‘after filtering’ data, I could explore the transfer function to see if the filtering is more involved than the description of the filter.

I will reply to you later when I have the opportunity to work on your data (and will delete this Comment).

[Markus Sause]

Here we go ...

File format is a little different as this is out of the commercial software. The signal was normalized before filtering(constant factor scaling), hence other amplitudes, first column is time in µs. Strong artefact at 150µs due to the step-function behaviour caused from zero-padding.

Many thanks for your effort!

[Star Strider]

Got it.

My pleasure!

[Star Strider]

Got it.

My pleasure!

--------------------------------------

Later ...

I cannot accurately identify that filter using the Signal Processing Toolbox or the System Identification Toolbox. A bandpass design comes close, but there is something the filter in your other software does that I cannot determine. Your filter takes the derivative of your signal, leading me to the bandpass design. Fine-tuning the bandpass filter will require experimentation.

The code I wrote to explore this is relatively straightforward. The ‘Ts’ variables are sampling interval, ‘Fs’ sampling frequency, ‘Fn’ Nyquist frequency, ‘Fv’ frequency vector, and ‘Iv’ index vector.

fidi1 = fopen('Markus Sause simwaves.txt','rt');
D1 = textscan(fidi1, '%f%f', 'CollectOutput',1);
Orig = cell2mat(D1);
Orig = Orig(1:end-1,:);
L1 = size(Orig,1);
fclose(fidi1);

fidi2 = fopen('Markus Sause simwave_filtered.txt','rt');        % Col #1 in µs
D2 = textscan(fidi2, '%f%f', 'CollectOutput',1);
Filt = cell2mat(D2);
L2 = size(Filt,1);
fclose(fidi2);

T1 = Orig(:,1);
Ts1 = mean(diff(T1));
Fs1 = 1/Ts1;
Fn1 = Fs1/2;
FT1 = fft(Orig(:,2))/L1;
Fv1 = linspace(0, 1, fix(L1/2)+1)*Fn1;
Iv1 = 1:length(Fv1);

T2 = Filt(:,1)*1E-6;
Ts2 = mean(diff(T2));
Fs2 = 1/Ts2;
Fn2 = Fs2/2;
FT2 = fft(Filt(:,2))/L2;
Fv2 = linspace(0, 1, fix(L2/2)+1)*Fn2;
Iv2 = 1:length(Fv2);

T3 = Orig(:,1)*1E-6;                                % Differentiated Signal
Ts3 = mean(diff(T2));
Fs3 = 1/Ts3;
Fn3 = Fs3/2;
Dorg = gradient(Orig(:,2),Ts1);
FT3 = fft(Dorg)/L2;
Fv3 = linspace(0, 1, fix(L1/2)+1)*Fn1;
Iv3 = 1:length(Fv2);

figure(1)
subplot(2,1,1)
plot(T1, Orig(:,2))
title('Original')
grid
subplot(2,1,2)
plot(T1, Filt(:,2))
set(gca, 'YLim',[-0.1  0.1])
title('Filtered')
grid

figure(2)
subplot(2,1,1)
plot(T1, Orig(:,2))
title('Original')
grid
subplot(2,1,2)
plot(T1, Dorg)
title('Derivative')
grid

TF23 = FT3./FT2;                             % Simple Transfer Function Approximation

figure(3)
subplot(3,1,1)
semilogy(Fv1, abs(FT1(Iv1))*2)
title('Original')
set(gca, 'XLim',[0 1E+6])
grid
subplot(3,1,2)
semilogy(Fv2, abs(FT2(Iv1))*2)
title('Filtered')
set(gca, 'XLim',[0 1E+6])
grid
subplot(3,1,3)
semilogy(Fv1, abs(TF23(Iv1))*2)
title('Transfer Function')
set(gca, 'XLim',[0 1E+6])
grid
xlabel('Frequency (Hz)')

Wp = [5E+4  4.5E+5]/Fn1;                                    % Specify Bandpass Filter
Ws = [5E+3  7.0E+5]/Fn1;
Rp = 10;
Rs = 30;
[n,Wn] = buttord(Wp, Ws, Rp, Rs);
[b,a] = butter(n,Wn);

[SOS,G] = tf2sos(b, a);                                     % Convert To SOS for Stability

figure(4)
freqz(b, a, 2^16, Fs1)

figure(5)
zplane(b, a)

DifFilt = filtfilt(SOS, G, Orig(:,2));

figure(6)
subplot(2,1,1)
plot(T1, Filt(:,2))
set(gca, 'YLim',[-0.1  0.1])
grid
subplot(2,1,2)
plot(T1, DifFilt)
grid

[Markus Sause]

Good morning!

First, many thanks for your efforts! You are reading my mind in assuming this software does something more, but it is hard to figure out.

Your SOS / filtfilt approach is exactly what I used to generate the very first images of this thread.

To me there are too many things I could imagine:

1. two-stage filter design (bandstop/bandpass)
2. subtraction of lines to obtain consistent periodicty
3. different choice of algorithm
4. Signal preprocessing (subsampling, padding)

Any recommendations, which are worthwhile to try?

[Star Strider]

My pleasure!

Give the filter an impulse and step response (separately) as an input and record the the output. That will reveal more about its characteristics than the narrow-band data you gave it. It then may be possible to identify it. Post those files here and I will see if it is possible to identify it.

Otherwise, without going through the code that does the filtering and reverse-engineering it (illegal because it is likely copyrighted or patented) there is no way to determine what the filtering algorithm actually is. I would write to the software author and ask them to describe in detail what the filtering algorithm does. Unless the filter is an important (rather than incidental) part of the software, there should not be a problem with this. It would then be possible to simulate it.

I have no idea what your signal is or what generates it. One option is to go the literature on it and see what is published on signal processing in whatever signal you are working with. (I do not know what you are doing, and I doubt that I have free access to those journals.) It is likely that the software implements an accepted and widely-used signal processing technique that was described in the literature before it was incorporated in the software.

[Markus Sause]

I really appreciate your help!

Attached is the filter response (same settings as before "6th order Butterworth, 100kHz) for a step-function 0->1 at digit 4000->4001 and for an impulse 0->1->0 with single digit zero at 4001. Both signals assume a 1e-7s sampling rate.

Of course I was asking the Software Company about the details. They told me many times they will let me know, but so far didn't. I do not think it is a company secret, but just lack of priority and organization. Regarding the science behind the signal, I am very well involved in the community and I know their standards. This is not typical and other software packages applied for this field are not capable of doing that. But no doubt, whatever has been implemented, it is likely been published someplace. But hard to search, if no name of the underlying procedure or algorithm. That is why I posted an image, hoping to find someone, who immediately recognizes the procedure.

[Star Strider]

My pleasure!

I’ll do what I can.

I downloaded them and will experiment with them later. These should allow us to identify it. It obviously involves a differentiator, since the step immediately falls back to baseline. It appears to be linear time invariant system.

[Star Strider]

I was able to identify the filter. There must be something else going on in that software, because implementing the filter with your original data does not reproduce the filtered data you want. There does not appear to be any delay, so there is no need to incorporate one.

I am surprised that the impulse and step response transfer functions are as different as they are. I expected some minor differences, but not the disparity I saw. Simulating the systems estimated from the impulse and step response data, with the impulse and step functions, demonstrated that the impulse response was not close to the original, while the step response was much closer. I would use the step response system ‘stpsys’ for any further studies you want to do with them.

The N (numerator) and D (denominator) transfer function coefficients the identification procedure estimates are:

Ni = [  0.000000000000000E+00    8.853420014635908E-01   -2.710878931309428E+00    1.883345694709287E+00    1.765616105760155E+00   -2.764546722715974E+00    9.411227595013788E-01  ];
Di = [  1.000000000000000E+00   -3.819282184119042E+00    5.475394472501968E+00   -3.494011282167925E+00    8.409925969220753E-01   -4.187786730629525E-03    1.124582487460075E-03  ];
Ns = [  0.000000000000000E+00   -5.186290325509300E-06    8.849533808620123E-01   -3.584355433029494E+00    5.444592985157094E+00   -3.675923276131520E+00    9.307375308562207E-01  ];
Ds = [  1.000000000000000E+00   -3.807722032895296E+00    5.453354497267004E+00   -3.509956946645446E+00    9.082260961276331E-01   -6.034062876276453E-02    1.648077149542787E-02  ];

The ‘i’ indicates the impulse response transfer function and the ‘s’ the step response transfer function.

The Code —

fidi1 = fopen('Markus Sause impulse-function_response.txt','rt');
imprespc = textscan(fidi1, '%f%f', 'CollectOutput',1);
impresp = cell2mat(imprespc);
impresp(:,1) = impresp(:,1)*1E-6;
fclose(fidi1);
fidi2 = fopen('Markus Sause step-function_response.txt','rt');
stprespc = textscan(fidi2, '%f%f', 'CollectOutput',1);
stpresp = cell2mat(stprespc);
stpresp(:,1) = stpresp(:,1)*1E-6;
fclose(fidi2);

figure(1)
subplot(2,1,1)
plot(impresp(:,1), impresp(:,2))
title('Impulse Response')
set(gca, 'XLim', [3.9  4.5]*1E-4)
grid
subplot(2,1,2)
plot(stpresp(:,1), stpresp(:,2))
title('Step Response')
set(gca, 'XLim', [3.9  4.5]*1E-4)
grid

Ts = mean(diff(impresp(:,1)));                                  % Sampling Interval
Fs = 1/Ts;                                                      % Sampling Frequency
Fn = Fs/2;                                                      % Nyquist Frequency
w = logspace(-1, 9, 1E+3);                                      % Frequency Vector (Hz)

zv = zeros(size(impresp(:,1)));                                 % Create Input Data
imp = zv;
imp(4000) = 1;
stp = zv;
stp(4000:end) = 1;

impdata = iddata(impresp(:,2), imp, Ts);                        % Create ‘iddata’ Objects
stpdata = iddata(stpresp(:,2), stp, Ts);

Np = 6;                                                         % Nr Poles
Nz = 6;                                                         % Nr Zeros

impsys = tfest(impdata, Np, Nz, 'Ts',Ts);                       % Estiamte Transfer Functions
Ni = impsys.Numerator;
Di = impsys.Denominator;
[magi,phasei,wouti,sdmagi,sdphasei] = bode(impsys, w);          % Calculate Bode Plot Data
stpsys = tfest(stpdata, Np, Nz, 'Ts',Ts);
Ns = stpsys.Numerator;
Ds = stpsys.Denominator;
[mags,phases,wouts,sdmags,sdphases] = bode(stpsys, w);

% fprintf(1, ['Ni = [' repmat('%23.15E  ', 1, Nz+1) '];\n'], Ni)
% fprintf(1, ['Di = [' repmat('%23.15E  ', 1, Np+1) '];\n'], Di)
% fprintf(1, ['Ns = [' repmat('%23.15E  ', 1, Nz+1) '];\n'], Ns)
% fprintf(1, ['Ds = [' repmat('%23.15E  ', 1, Np+1) '];\n'], Ds)

figure(2)
subplot(2,1,1)
semilogx(wouti, 20*log10(squeeze(magi)))
title('Impulse Response Bode Plot')
ylabel('Magnitude (dB)')
axis([min(w)  max(w)    -100  50])
grid
subplot(2,1,2)
semilogx(wouti, squeeze(phasei))
xlabel('Frequency (Hz)')
ylabel('Phase (°)')
axis([min(w)  max(w)    0  1200])
grid

figure(3)
subplot(2,1,1)
semilogx(wouts, 20*log10(squeeze(mags)))
title('Step Response Bode Plot')
ylabel('Magnitude (dB)')
axis([min(w)  max(w)    -100  50])
grid
subplot(2,1,2)
semilogx(wouts, squeeze(phases))
xlabel('Frequency (Hz)')
ylabel('Phase (°)')
axis([min(w)  max(w)    0  1200])
grid

fidi3 = fopen('Markus Sause simwaves.txt','rt');                % Load Original Data
D1 = textscan(fidi1, '%f%f', 'CollectOutput',1);
Orig = cell2mat(D1);
Orig = Orig(1:end-1,:);
L1 = size(Orig,1);
fclose(fidi3);

[sos,g] = tf2sos(Ns, Ds);                                       % Convert To Second Order Section

figure(4)
freqz(sos, 2^16, Fs)                                            % Analyse Filter
set(subplot(2,1,1), 'XScale','log')
set(subplot(2,1,2), 'XScale','log')

OrigFilt = filtfilt(sos, g, Orig(:,2));                         % Filter & Plot

figure(5)
plot(Orig(:,1), OrigFilt)
grid

figure(6)
impulse(impsys)

figure(7)
step(stpsys)

This assumes you have the System Identification Toolbox. If you don’t, and want me to post the plots, I will.

[Markus Sause]

Sorry for being late on this!

Many thanks! I will give it a try and let you know my feedback asap.

[Star Strider]

My pleasure!

No worries! We all have lives outside MATLAB Answers!

[Star Strider]

My pleasure!

Consider using filtfilt instead of filter. The filtfilt function has a maximally-flat phase response, so it will not produce any phase distortion in your filtered signal.

If my Answer helped you solve your problem, please Accept it!

[Markus Sause]

Surprisingly the filtfilt destroys the performance of the filter, so I used the filter command instead.

Thanks!

[Star Strider]

My pleasure!

The filtfilt behavior surprises me as well.

Answer 2

[Markus Sause]

Star Strider,

you showed me the way to go! Based on your code I tried and evaluated several filter options. With the estimates from impulse and step function I concluded for some of the parameters.

What comes closest is:

filterfreq = 100e3;  
hpFilt = designfilt('highpassiir','StopbandFrequency',filterfreq/4, ...
         'PassbandFrequency',filterfreq,...
         'StopbandAttenuation',80,'DesignMethod','butter', 'SampleRate', Fs);
fvtool(hpFilt);
OrigFilt = filter(hpFilt, Orig(:,2));

Many thanks for your help!

Markus