Which filtering approach to use?
Show older comments
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
2 Comments
Chad Greene
on 27 Feb 2017
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
on 28 Feb 2017
Accepted Answer
Star Strider
on 28 Feb 2017
1 vote
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.
15 Comments
Markus Sause
on 28 Feb 2017
Edited: Markus Sause
on 28 Feb 2017
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
on 28 Feb 2017
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
on 28 Feb 2017
Star Strider
on 28 Feb 2017
Got it.
My pleasure!
Star Strider
on 1 Mar 2017
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
on 1 Mar 2017
Star Strider
on 1 Mar 2017
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
on 2 Mar 2017
Star Strider
on 2 Mar 2017
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
on 3 Mar 2017
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
on 6 Mar 2017
Star Strider
on 6 Mar 2017
My pleasure!
No worries! We all have lives outside MATLAB Answers!
Star Strider
on 14 Mar 2017
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
on 15 Mar 2017
Star Strider
on 15 Mar 2017
My pleasure!
The filtfilt behavior surprises me as well.
More Answers (1)
Markus Sause
on 14 Mar 2017
Categories
Find more on Digital Filter Design in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!
