Gibbs-like behaviour with lowpass on long signals

I have a signal that is 1000 some data points sampled at a rate of 1/3 hz. I'd like to lowpass it with fpass = 0.08.

I have an example attached the first 196 of them:

t=readtable('example_signal.txt');
filtered=lowpass(t.Var1,0.08,1/3); 
plot(filtered)

But this gives something that looks like gibbs!

t=readtable('example_signal.txt');
filtered=lowpass(t.Var1([1:178]),0.08,1/3); 
plot(filtered)

And this doesn't?

Why does making the signal shorter change things? If I increase steepness, I can avoid this, but if the signal is longer yet again, I need to keep increasing steepness, which I can only do so much of -- and I can't make it steep enough to do the whole 1k long signal without an issue.

filtered=lowpass(t.Var1,0.08,1/3,'Steepness',0.9); 
plot(filtered)

My guess is that as the signal lengthens, some part of the Fourier transform of it is getting pushed into the transition band, and that's why making it steeper helps.

But I don't know why it's happening, or how to avoid it.

If I knew how to calculate when this would happen, then I could just break my signal into parts shorter than this critical threshold. But I don't understand why that should help theoretically (indeed, it seems like I'd be throwing out information on very slow oscillations (which should pass anyways) if I were to do that).

Answer 1

[Star Strider]

Your code and results appear to be appropriate. The only change I made here was to add 'ImpulseResponse='iir''

to the lowpass calls. The lowpass (and other such functions) automatically design FIR filters (in most instances). Overriding this by specify in IIR design (an efficient elliptic filter) usually gives better results. Also, subtracting the D-C offset from the original signal sometimes gives better results, because of the significant trtansient the D-C offset would otherwise cause. (That did not appear to make a difference here.)

Try something like this --

t = readtable('example_signal.txt');
sv = t.Var1;                                                % Signal Vector
Fs = 1/3;                                                   % Sampling Frequency (Hz)
tv = linspace(0, numel(sv)-1, numel(sv))/Fs;                % Time Vector

svm = mean(sv);

figure
plot(tv, sv)
grid
xlabel('Time (s)')
ylabel('Amplitude (units)')
title('Original Signal')

Fco = 0.08;

[FTs1,Fv] = FFT1(sv,tv);

figure
plot(Fv, abs(FTs1)*2)
grid
xlabel('Frequency (Hz)')
ylabel('Magnitude (units)')
title('Fourier Transform Of Original Signal')
xline(Fco, '--r', 'F_{co}')

filtered=lowpass(t.Var1,0.08,1/3, ImpulseResponse='iir'); 

figure
plot(tv, filtered)
grid
xlabel('Time (s)')
ylabel('Amplitude (units)')
title('Filtered Signal (D-C Offset Retained)')

●

filtered2 = lowpass(t.Var1-svm,0.08,1/3, ImpulseResponse='iir'); 

figure
plot(tv, filtered2+svm)
grid
xlabel('Time (s)')
ylabel('Amplitude (units)')
title('Filtered Signal (D-C Offset Subtracted, Then Added Back)')

[FTs1,Fv] = FFT1(filtered,tv);

figure
plot(Fv, abs(FTs1)*2)
grid
xlabel('Frequency (Hz)')
ylabel('Magnitude (units)')
title('Fourier Transform Of Filtered Signal')
xline(Fco, '--r', 'F_{co}')

function [FTs1,Fv] = FFT1(s,t)
% One-Sided Numerical Fourier Transform
% Arguments:
%    s: Signal Vector Or Matrix
%    t: Associated Time Vector 
    t = t(:);
    L = numel(t);
    if size(s,2) == L
        s = s.';
    end
    Fs = 1/mean(diff(t));
    Fn = Fs/2;
    NFFT = 2^nextpow2(L);
    FTs = fft((s - mean(s)) .* hann(L).*ones(1,size(s,2)), NFFT)/sum(hann(L));
    Fv = Fs*(0:(NFFT/2))/NFFT;
%     Fv = linspace(0, 1, NFFT/2+1)*Fn;
    Iv = 1:numel(Fv);
    Fv = Fv(:);
    FTs1 = FTs(Iv,:);
end

Calculating the Fourier transform also usually helpt to clarify these sorts of problems.

.

Answer 2

[Paul]

Hi qfn,

I think that lowpass, et. al. will result in transients due the large, effective steps at the endpoints.

I think this issue (or something like it) is dicussed in the documentation or in other threads here on Answers, though I can't find either at the moment.

In any event, for this particular data set, subtracting the mean -> lowpass -> add the mean seems to work pretty well:

T=readtable('example_signal.txt');
x = T.Var1;
Fs = 1/3;
t = (0:numel(x)-1)/Fs;
[xf1,d] = lowpass(x,0.08,Fs);  % original
xf2 = lowpass(x-mean(x),0.08,Fs) + mean(x);
figure
plot(t,x,t,xf1*nan,t,xf2),grid  % removed xf1 for clarity

●