How to remove the ripple or artifacts caused by lowpass filter?

Hi, I am in a complex situation, where the lowpass fitering with sharp cutoff freqeuncy is affecting a technique I was developing. The following first code use a fc = 80Hz, to extract the envelop. As shown in the second code, I want to reduce the fc to 15Hz and achieve a smooth envelop.i.e. gradient A' reduces in the second simulation compared to first due to reduction of fc from 80 Hz to 15 Hz.

But, reducing to 15Hz, creates additional ripples and sharp edges and changes the expected envelop, which is a problem for me. So,Is there any other way to escape from these artifacts? I tried reducing the filter order to 1, but butter filtering is not performing well. So I am looking if there is any other way to elminate those added artifacts either by changing the lowpass filteration or by post proceesing the extracted envelop?

clear all; 
% close all; clc;
%-------------------------------------------------------------
% System definition %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%-------------------------------------------------------------

fs=8192;        % sampling frequency                                       
dt = 1/fs;       % sample time                                             
T=3;             % duration of the signal                                       
Nt = T*fs;       % total number of samples
t = 0:dt:T-dt;   % time vector

% Source definition
f0 = 900;    % initial angular speed in Hz
fT = 2700;    % final angular speed in Hz
finst = linspace(f0,fT,Nt);
phi = 2*pi*cumsum(finst)*dt;


a1 = 45; 
b1 = 5000; 
c1 = 2; 
c2 = 1.5;
d1 = 2+3*1j;

% Definition of envelop
int_phi0 = 3*pi/4;
A = d1+(a1 * exp(-b1 * (t - T/c1).^2 + 1j * int_phi0));


% Definition of source signal
q = A.'.*exp(1j*phi).';


filtorder =6; fc =80;%----------------------------------% fc-cutoff
[bcoeff,acoeff] = butter(filtorder,fc/fs*2);
Afilt = filtfilt(bcoeff,acoeff,q.*exp(-1j*phi).');


% figure()
% freqz(bcoeff,acoeff,[],fs)

figure()
plot(t,q,'k')

Warning: Imaginary parts of complex X and/or Y arguments ignored.

hold on
plot(t,abs(Afilt),'r','LineWidth',2)
grid on; box on;
legend('q(t)','A(t)')
xlabel('s')
ylabel('Amp')
title('q(t) and A(t) - filtered')
set(gca,'FontSize',15)
xlim([1.2 1.8])

subplot()
clear all; 
% close all; clc;
%-------------------------------------------------------------
% System definition %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%-------------------------------------------------------------

fs=8192;        % sampling frequency                                       
dt = 1/fs;       % sample time                                             
T=3;             % duration of the signal                                       
Nt = T*fs;       % total number of samples
t = 0:dt:T-dt;   % time vector

% Source definition
f0 = 900;    % initial angular speed in Hz
fT = 2700;    % final angular speed in Hz
finst = linspace(f0,fT,Nt);
phi = 2*pi*cumsum(finst)*dt;


a1 = 45; 
b1 = 5000; 
c1 = 2; 
c2 = 1.5;
d1 = 2+3*1j;

% Definition of envelop
int_phi0 = 3*pi/4;
A = d1+(a1 * exp(-b1 * (t - T/c1).^2 + 1j * int_phi0));


% Definition of source signal
q = A.'.*exp(1j*phi).';


filtorder =6; fc =15;%----------------------------------% fc-cutoff
[bcoeff,acoeff] = butter(filtorder,fc/fs*2);
Afilt = filtfilt(bcoeff,acoeff,q.*exp(-1j*phi).');


% figure()
% freqz(bcoeff,acoeff,[],fs)

figure()
plot(t,q,'k')

Warning: Imaginary parts of complex X and/or Y arguments ignored.

hold on
plot(t,abs(Afilt),'r','LineWidth',2)
grid on; box on;
legend('q(t)','A(t)')
xlabel('s')
ylabel('Amp')
title('q(t) and A(t) - filtered')
set(gca,'FontSize',15)
xlim([1.2 1.8])

Answer 1

[Mathieu NOE]

hello

Butterworth filters do have ripples in their step response...

look at the step response of your Butterworth filter with fc = 80 Hz :

(x axis is samples not seconds as displayed !)

now that you have lowered the fc to 15 Hz, it's getting very long until the filter's output is steady again :

no surprise , it takes 5 time longer and with oscillations

so you should maybe use Bessel filters instead (step response do not exhibit any oscillation)

this is my suggestion below , you can change filtorder and fc as you wish, this will affect how close the envelop will be to the "real" one, but the ripples are gone

% clear all; 
% close all; clc;
%-------------------------------------------------------------
% System definition %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%-------------------------------------------------------------
fs=8192;        % sampling frequency                                       
dt = 1/fs;       % sample time                                             
T=3;             % duration of the signal                                       
Nt = T*fs;       % total number of samples
t = 0:dt:T-dt;   % time vector
% Source definition
f0 = 900;    % initial angular speed in Hz
fT = 2700;    % final angular speed in Hz
finst = linspace(f0,fT,Nt);
phi = 2*pi*cumsum(finst)*dt;
a1 = 45; 
b1 = 5000; 
c1 = 2; 
c2 = 1.5;
d1 = 2+3*1j;
% Definition of envelop
int_phi0 = 3*pi/4;
A = d1+(a1 * exp(-b1 * (t - T/c1).^2 + 1j * int_phi0));
% Definition of source signal
q = A.'.*exp(1j*phi).';
% filtorder =6; fc =80;%----------------------------------% fc-cutoff
% [bcoeff,acoeff] = butter(filtorder,fc/fs*2);

% alternative with Bessel filter 
filtorder =3; fc =500;%----------------------------------% fc-cutoff
[b,a] = besself(filtorder,fc,'low'); % NB it's analog ! 
[bcoeff,acoeff] = c2dm(b,a,1/fs,'tustin');

Afilt = filtfilt(bcoeff,acoeff,q.*exp(-1j*phi).');
% figure()
% freqz(bcoeff,acoeff,[],fs)
figure()
dstep(bcoeff,acoeff)

figure()
plot(t,q,'k')
hold on
plot(t,abs(Afilt),'r','LineWidth',2)
grid on; box on;
legend('q(t)','A(t)')
xlabel('s')
ylabel('Amp')
title('q(t) and A(t) - filtered')
set(gca,'FontSize',15)
xlim([1.2 1.8])
% 
% subplot()

clear all; 
% close all; clc;
%-------------------------------------------------------------
% System definition %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%-------------------------------------------------------------
fs=8192;        % sampling frequency                                       
dt = 1/fs;       % sample time                                             
T=3;             % duration of the signal                                       
Nt = T*fs;       % total number of samples
t = 0:dt:T-dt;   % time vector
% Source definition
f0 = 900;    % initial angular speed in Hz
fT = 2700;    % final angular speed in Hz
finst = linspace(f0,fT,Nt);
phi = 2*pi*cumsum(finst)*dt;
a1 = 45; 
b1 = 5000; 
c1 = 2; 
c2 = 1.5;
d1 = 2+3*1j;
% Definition of envelop
int_phi0 = 3*pi/4;
A = d1+(a1 * exp(-b1 * (t - T/c1).^2 + 1j * int_phi0));
% Definition of source signal
q = A.'.*exp(1j*phi).';
% filtorder =6; fc =15;%----------------------------------% fc-cutoff
% [bcoeff,acoeff] = butter(filtorder,fc/fs*2);

% alternative with Bessel filter 
filtorder =3; fc =100;%----------------------------------% fc-cutoff
[b,a] = besself(filtorder,fc,'low'); % NB it's analog ! 
[bcoeff,acoeff] = c2dm(b,a,1/fs,'tustin');

Afilt = filtfilt(bcoeff,acoeff,q.*exp(-1j*phi).');

% figure()
% freqz(bcoeff,acoeff,[],fs)
figure()
dstep(bcoeff,acoeff)

figure()
plot(t,q,'k')
hold on
plot(t,abs(Afilt),'r','LineWidth',2)
grid on; box on;
legend('q(t)','A(t)')
xlabel('s')
ylabel('Amp')
title('q(t) and A(t) - filtered')
set(gca,'FontSize',15)
xlim([1.2 1.8])

[Kalasagarreddi Kottakota]

Hi @Mathieu NOE,

Thanks, Thats indeed very helpfull. But i have question, where 'fc' is really cuttoff frequency of lowpass filter?

As, In matlab, it is defined Wo as the frequency Wo is the frequency up to which the group delay is approximately constant.

[B,A] = besself(N,Wo)

[Mathieu NOE]

oops sorry

I passed the fc in Hz whereas it should be in rad /s

change the lines

[b,a] = besself(filtorder,fc,'low'); % NB it's analog !

into

[b,a] = besself(filtorder,2*pi*fc,'low'); % NB it's analog !

and divide fc values I used by 2pi to keep output the same

[Kalasagarreddi Kottakota]

Hi @Mathieu NOE, I made the change of 2pi, but when I plot the freqeuncy response, it is very unpredictable, fc = 500 in the below code, what can be read from this freqeuncy response?

fs = 8192;
filtorder =3; fc =500;%----------------------------------% fc-cutoff
[b0,a0] = besself(filtorder,2*pi*fc,'low'); % NB it's analog ! 
[bcoeff,acoeff] = c2dm(b0,a0,1/fs,'tustin');

figure()
freqz(bcoeff,acoeff,[],fs)
xlim([0 1000])

[Star Strider]

Bessel filters cannot be represented as digital filters, and retain their phase-neutral characteristics.

[Mathieu NOE]

@Kalasagarreddi Kottakota

I don't see what is the trouble with the Bode plot of the bessel filter

(also done now in your code below)

at f = fc the filter has a -6 dB attenuation

also this is your code , now corrected . notice that with the 2pi corrction we finally are using the same fc as in your original code . It's just now bessel instead of butter and the filter order has been reduced from 6 to 3

clear all; 
close all; clc;
%-------------------------------------------------------------
% System definition %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%-------------------------------------------------------------
fs=8192;        % sampling frequency                                       
dt = 1/fs;       % sample time                                             
T=3;             % duration of the signal                                       
Nt = T*fs;       % total number of samples
t = 0:dt:T-dt;   % time vector
% Source definition
f0 = 900;    % initial angular speed in Hz
fT = 2700;    % final angular speed in Hz
finst = linspace(f0,fT,Nt);
phi = 2*pi*cumsum(finst)*dt;
a1 = 45; 
b1 = 5000; 
c1 = 2; 
c2 = 1.5;
d1 = 2+3*1j;
% Definition of envelop
int_phi0 = 3*pi/4;
A = d1+(a1 * exp(-b1 * (t - T/c1).^2 + 1j * int_phi0));
% Definition of source signal
q = A.'.*exp(1j*phi).';
% filtorder =6; fc =80;%----------------------------------% fc-cutoff
% [bcoeff,acoeff] = butter(filtorder,fc/fs*2);
% alternative with Bessel filter 
filtorder =3; fc =80;%----------------------------------% fc-cutoff
[b,a] = besself(filtorder,2*pi*fc,'low'); % NB it's analog ! 
[bcoeff,acoeff] = c2dm(b,a,1/fs,'tustin');
Afilt = filtfilt(bcoeff,acoeff,q.*exp(-1j*phi).');
figure()
f = logspace(1,3,300);
freqz(bcoeff,acoeff,f,fs)
figure()
dstep(bcoeff,acoeff)
figure()
plot(t,q,'k')
hold on
plot(t,abs(Afilt),'r','LineWidth',2)
grid on; box on;
legend('q(t)','A(t)')
xlabel('s')
ylabel('Amp')
title('q(t) and A(t) - filtered')
set(gca,'FontSize',15)
xlim([1.2 1.8])
% 
% subplot()
clear all; 
% close all; clc;
%-------------------------------------------------------------
% System definition %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%-------------------------------------------------------------
fs=8192;        % sampling frequency                                       
dt = 1/fs;       % sample time                                             
T=3;             % duration of the signal                                       
Nt = T*fs;       % total number of samples
t = 0:dt:T-dt;   % time vector
% Source definition
f0 = 900;    % initial angular speed in Hz
fT = 2700;    % final angular speed in Hz
finst = linspace(f0,fT,Nt);
phi = 2*pi*cumsum(finst)*dt;
a1 = 45; 
b1 = 5000; 
c1 = 2; 
c2 = 1.5;
d1 = 2+3*1j;
% Definition of envelop
int_phi0 = 3*pi/4;
A = d1+(a1 * exp(-b1 * (t - T/c1).^2 + 1j * int_phi0));
% Definition of source signal
q = A.'.*exp(1j*phi).';
% filtorder =6; fc =15;%----------------------------------% fc-cutoff
% [bcoeff,acoeff] = butter(filtorder,fc/fs*2);
% alternative with Bessel filter 
filtorder =3; fc =15;%----------------------------------% fc-cutoff
[b,a] = besself(filtorder,2*pi*fc,'low'); % NB it's analog ! 
[bcoeff,acoeff] = c2dm(b,a,1/fs,'tustin');
Afilt = filtfilt(bcoeff,acoeff,q.*exp(-1j*phi).');
figure()
f = logspace(0,2,100);
freqz(bcoeff,acoeff,f,fs)
figure()
dstep(bcoeff,acoeff)
figure()
plot(t,q,'k')
hold on
plot(t,abs(Afilt),'r','LineWidth',2)
grid on; box on;
legend('q(t)','A(t)')
xlabel('s')
ylabel('Amp')
title('q(t) and A(t) - filtered')
set(gca,'FontSize',15)
xlim([1.2 1.8])

Answer 2

[Walter Roberson]

Imagine that you were doing a discrete fourier transform of a square wave. Maybe it looks not bad at a target frequency. But then you try to use it at a higher frequency, and you look carefully at the leading and trailing edges and see ripples. Would decreasing the number of fourier coefficients solve the problem? Obviously not. Would increasing the number of fourier coefficients solve the problem? No -- examine the result at a sufficiently high frequency and you would continue to see ripples at the leading and trailing edges.

Using butter filter has similar problems; it is a discrete reconstruction of an ideal filter shape, and if you use too few coefficients the result is going to look bad. You can increase the filter order to make the ripples less obvious, but you cannot eliminate them with any finite order. The best you can do is reduce the magnitude of the ripple to less than a level that is acceptable to you -- which might require increasing the order (which in turn implies increased delays). buttord might help in the design.