Low pass filtering in software

[gf]

As I understand, it happens due to Gibbs phenomenon because filter removes some components from signal. But how to estimate its worse case?

The Gibbs Phenomenon basically says that a square wave cannot be represented exactly by its Fourier series, not even with an an infinite number of harmonics. With increasing number of included harmonics, the area between the approximation and the true square wave does indeed converge to zero, but the maximum absolute deviation does not. An overshoot always remains - it just becomes narrower and narrower on the time axis.

There is a closed-form solution for the asymptotic limit of the overshoot amplitude, and it is the Wilbraham-Gibbs Constant: 1.17897974447216 ...

[radiolistener]

If you want the maximum response from any input consisting only of ±1, construct one that matches the signs of the input filter in reverse order.  Then, the maximum amplitude is the sum of the magnitudes of the filter coefficients,
$$\sum_{k=0}^{N-1} \lvert c_k \rvert$$
and occurs when the input signs match the corresponding coefficient signs.

Cool, it works

sum(abs(kernel)): 2.34858



But 2.34858 overshoot amplitude looks too much. Is there a way to minimize it during filter design?

[gf]

Here is a small test code in Octave/Matlab which shows issue with signal overshoot on the FIR output:

Code: [Select]

pkg load signal;

Fs = 48000*32;     % Sample rate [Hz]
Fpass = 22000;     % Passband [Hz]

Wp = Fpass / (Fs / 2);  % Normalized passband
beta = 6.181877;        % Kaiser window beta (60 dB stopband attenuation, 0.01 dB passband ripple)
N = 3040;               % FIR order (estimated with ceil((A - 8) / (2.285 * pi * Wt)))

kernel = fir1(N, Wp, kaiser(N + 1, beta));


I assume that you want the passband extend up to 22kHz (within the desired passband ripple).
However, the fir1 parameter Wn is not the end of the passband, but approx. the center of the transition band.
I suggest to use kaiserord() to calculate Wn, beta and N, e.g.

Code: [Select]

passband_ripple = 0.001  % 20*log10(1+0.001) = 0.00868 dB
stopband_atten = 0.001   % 60 dB
Fpass = 22000    % end of passband
Fstop = 24000    % begin of stopband
[n,Wn,beta,ftype] = kaiserord([Fpass Fstop], [1 0],  [passband_ripple stopband_atten], Fs)
kernel = fir1(n,Wn,ftype,kaiser(n+1,beta));


[radiolistener]

I assume that you want the passband extend up to 22kHz (within the desired passband ripple).
However, the fir1 parameter Wn is not the end of the passband, but approx. the center of the transition band.
I suggest to use kaiserord() to calculate Wn, beta and N, e.g.

I simplified parameters selection for this test code, so it don't include its estimation.

Regarding the kaiserord function, I know about it, I was started to design FIR with Kaiser window with using it. But it appears that it is very problematic. I don't remember all details, but I even wanted to open Octave issue about mistakes in this function. As I remember it has many issues, include different results with Matlab and its code don't corresponds with classic filter design using Kaiser window explained in science articles. As I remember there is some mistakes in its code (or may be improvements, but they are missing in original article). So I'm avoid to use it.

[Tation]

No matter what "classical" method you use to design this FIR filter (windowed with any window, Parks-McLellan,...) it will show overshot because of the Gibbs phenomenon.

But there are other design methods that allow the specification of both frequency and tine domain constraints, being impulse overshot a common constraint. Remember seeing some IEEE papers on that some years ago, but do not know if Octave/Matlab or any other software tool for digital filter design implement any of this methods.

[tggzzz]

In order to preserve waveform shape, obviously a critical point for an oscilloscope, decent analogue scopes' front ends have a Gaussian response.

[SiliconWizard]

No matter what "classical" method you use to design this FIR filter (windowed with any window, Parks-McLellan,...) it will show overshot because of the Gibbs phenomenon.

But there are other design methods that allow the specification of both frequency and tine domain constraints, being impulse overshot a common constraint. Remember seeing some IEEE papers on that some years ago, but do not know if Octave/Matlab or any other software tool for digital filter design implement any of this methods.

For instance, averaging will not overshoot by nature, and a simple moving average is actually a particular type of FIR where all coefficients are equal, positive and with a 1/N value (with N the number of coefficients = size of the moving average window). A moving average has characteristics that may not fit a particular use case though, like relatively large side-lobes past Fs / N (where Fs = sampling frequency).

[Psi]

And if you want to get fancy, have a moving average that dynamically changes its depth based on a rate of change calculation.
I did that once for some accelerometer stuff and it worked amazingly well. Not saying it was the best choice, but it was easy to understand and didn't need any complex/slow math.  It responded fast to big changes but also smoothed out idle noise very well

[SiliconWizard]

For preserving signal extrema while still low-pass filtering, I've otherwise used Savitzky-Golay filters. (Note that the moving average is a simple form of Savitzky-Golay filter.)

[fourfathom]

Yes, a FIR does not necessarily have overshoot.  I needed to implement a Gaussian-response FIR for a frequency shift keying (8GFSK) design and someone here pointed me to a nice implementation.  This one doesn't use the classical FIR with coefficients design, but instead a chained array of 2-stage boxcar filters.  Here's a link to my description of it: [url]http://wb6cxc.com/?p=126]http://wb6cxc.com/?p=126] [url]http://wb6cxc.com/?p=126[/url]

Anyway, it's a FIR without ringing.  It uses more delay registers than a classical implementation, but it doesn't need the coefficient storage and multiply stages, and lets me easily adjust the filter length when changing GFSK modes.

[gf]

Regarding the kaiserord function, I know about it, I was started to design FIR with Kaiser window with using it. But it appears that it is very problematic. I don't remember all details, but I even wanted to open Octave issue about mistakes in this function. As I remember it has many issues, include different results with Matlab and its code don't corresponds with classic filter design using Kaiser window explained in science articles. As I remember there is some mistakes in its code (or may be improvements, but they are missing in original article). So I'm avoid to use it.

Basically it appears to work, but to be honest I have not tested it thoroughly. Based on your experience, it seems I should use it with care, not trusting the results blindly.

[benseno]

I saw that Gibbs Phenomenon was mentioned. Isn't it true only for cases when signal is converted to the frequency domain and back to the time domain? FIR doesn't do that indeed.

[Nominal Animal]

No. Consider a pulse train of 2π wavelength, first half of each wave at +1 and latter half at -1.  Any periodic signal can be represented by an infinite series
$$y(t) = a_0 + \sum_{n=0}^{\infty} \biggl( a_n \cos( n t ) + b_n \sin( n t ) \biggr)$$
where
$$\left\lbrace ~ \begin{aligned}
a_0 &= \frac{1}{P} \int_{0}^{P} y(t) \, d t \\
a_n &= \frac{2}{P} \int_{0}^{P} y(t) \cos (n t) \, d t, \quad n \ge 1 \\
b_n &= \frac{2}{P} \int_{0}^{P} y(t) \sin (n t) \, d t, \quad n \ge 1 \\
\end{aligned} \right .$$
where \$P\$ is the signal period.  If you work out the above, you find that for this signal \$P = 2 \pi\$ and
$$y(t) = \frac{4}{\pi} \sum_{k=0}^{\infty} \frac{1}{2 k + 1} \sin\bigl((2 k + 1) t\bigr)$$
Essentially, it can be represented as a sine series with only odd \$n\$ of the base frequency have a nonzero amplitude, \$b_n = 4/(\pi n)\$, with all others zero (\$a_0 = a_n = 0\$):
$$y(t) = \frac{4}{\pi} \sin(t) + \frac{4}{3 \pi} \sin(2 t) + \frac{4}{5 \pi} \sin(3 t) + \frac{4}{7 \pi} \sin(5 t) + \dots$$
In simple terms, each higher frequency component cancels out most of the overshoot of the previous frequencies, moving the overshoot closer to the edge.  (In case you wonder, with \$N\$ first terms of the sine series, the peaks occur when \$\sum_{n=0}^{N-1} \cos \bigl((2 n + 1) x\bigr) = 0\$.)

This means that any filter with a rectangular frequency response will suffer from Gibbs' Phenomenon.

This is because you just chop off the equivalent series representation at some \$n\$, when you use a rectangular frequency response filter.  As shown above, that will result in overshoot (Gibbs' Phenomenon) at maximal slew rate edges.  This is physically real, too: you see this in practical analog filters that only remove highest frequencies and do not amplify any frequencies, as well.

(Moreover, as the amplitude of the \$n\$'th term of the example pulse train's Fourier transform is \$4/((2 n + 1) \pi)\$, it means the amplitude of the component sine waves drop linearly: if you want to double the precision of your approximation, you need double the number of sine waves!  The convergence is annoyingly slow.  For real life, it means you need a surprisingly wide spectrum response to reproduce rectangular pulse trains at high fidelity.  Or, more simply, that those pesky high-slewrate edges have surprisingly high frequency components in them.)

To avoid Gibbs' Phenomenon, you cannot have a rectangular frequency response: you need to attenuate the highest frequencies.

[gf]

Strictly speaking, the Gibbs Phenomenon is concerned with the under/overshoot near the discontinuity points (let me call them "Gibbs points") when a discontinuous periodic signal is resynthesised from its Fourier series. It is nicely explained e.g. here.

Regarding FIR filter: The samples of a discrete-time signal always represent a bandwidth-limited signal in the continuous-time domain, and a bandwidth-limited signal is never discontinuous. So strictly, I don't see how a discrete-time signal can have any Gibbs points.

Nevertheless, similar under/overshoots can, of course, also be observed near smoothed (continuous) edges as long as they are steep enough and if the signal is resynthesised/approximated with truncated Fourier series. But I find it difficult to decide whether we can still call this the Gibbs Phenomenon or not. It's a question of how strictly you interpret the definition. The difference between such edges an true Gibbs points is that the under/overshoot amplitude near smooth edges converges to zero as we include more and more terms of the Fourier series. At true Gibbs points (i.e. discontinuities) the overshoot amplitude does not converge to zero even if we include an infinite number of terms.

[radiolistener]

Your point about a discrete-time signal always representing a bandwidth-limited continuous-time signal is correct. However, it's important to consider the role of an FIR filter in this context.

Although an FIR filter operates on discrete samples, it effectively performs convolution with a discrete-time impulse response, which can be interpreted as processing a continuous-time signal. This means that the output samples of the FIR filter no longer correspond to the aperture of the original discrete-time signal but rather to a newly sampled version of the underlying continuous-time representation.

Thus, while a properly band-limited signal does not inherently exhibit Gibbs points, FIR filtering can introduce effects that resemble them, particularly when dealing with signals with sharp edges or when using truncated filter kernels.

[SiliconWizard]

Regarding FIR filter: The samples of a discrete-time signal always represent a bandwidth-limited signal in the continuous-time domain, and a bandwidth-limited signal is never discontinuous. So strictly, I don't see how a discrete-time signal can have any Gibbs points.

You may be missing one point when it comes to FIR depending on how you design them. The impulse response truncation method is bound to exhibit the Gibbs phenomenon:

https://www.gaussianwaves.com/2010/04/gibbs-phenomena-a-demonstration/

[Nominal Animal]

At true Gibbs points (i.e. discontinuities) the overshoot amplitude does not converge to zero even if we include an infinite number of terms.

Yes; it is just the overshoot width or duration, not amplitude, that tends to zero.

However, we don't have a better name for this Gibbs-like phenomenon; it's not exactly ringing or plain overshoot either.

[/quote]

The impulse response truncation method is bound to exhibit the Gibbs phenomenon:

  This is exactly what I tried to say in #37, using way too many words.

[benseno]

Will the Gibbs phenomenon be appeared in the result of backward Fourier transform of a discontinuous signal (in the frequency domain) after frequency domain equivalent of smoothing window filter is applied to it?

[Nominal Animal]

Time domain and frequency domain are just two different ways of describing the same thing, they're not separate things.

If you take any full-spectrum signal with sharp edges, then just remove all the highest frequencies, you will see the same kind of overshoot as in Gibbs phenomenon.  This applies to both the math, and to physical signals.

[benseno]

Seems i couldn't formulate my question properly (thanks to my bad English).
Let assume 2 methods of filtering:
(a) directly by FIR filter
(b) forward Fourier transform is done to discontinuous square wave signal, then equivalent of FIR filter in the frequency domain is applied, and last step - backward Fourier transform is done to rebuilt signal on the time domain.

The question is - will the Gibbs phenomenon occur in method (b) if (a) is not resulted in it?

Thanks

[SiliconWizard]

Seems i couldn't formulate my question properly (thanks to my bad English).
Let assume 2 methods of filtering:
(a) directly by FIR filter
(b) forward Fourier transform is done to discontinuous square wave signal, then equivalent of FIR filter in the frequency domain is applied, and last step - backward Fourier transform is done to rebuilt signal on the time domain.

The question is - will the Gibbs phenomenon occur in method (b) if (a) is not resulted in it?

Thanks

Method (b) will give you similar issues if you use a rectangular window.
To simplify the concept, in the case of a FIR filter (designed with the impulse response truncation), the discontinuity can be found in the step function.
Using convolution (direct Fourier transform, weights in the frequency domain, inverse Fourier transform), the discontinuity will be found in the window if it's rectangular.

So method (b) requires using appropriate windows.
Note that unless you need to filter only a relatively short buffer (and you don't need to do it in real time), using a single direct and inverse FT on the whole signal becomes unpractical and the usual way of doing it is the overlap-add method (and in particular the WOLA), in which you process a sequence of overlapping sections of signal:

https://www.dsprelated.com/freebooks/sasp/Overlap_Add_OLA_STFT_Processing.html

 

[gf]

Seems i couldn't formulate my question properly (thanks to my bad English).
Let assume 2 methods of filtering:
(a) directly by FIR filter
(b) forward Fourier transform is done to discontinuous square wave signal, then equivalent of FIR filter in the frequency domain is applied, and last step - backward Fourier transform is done to rebuilt signal on the time domain.

The question is - will the Gibbs phenomenon occur in method (b) if (a) is not resulted in it?

No, if the resulting signal (after convolution with the FIR kernel) has no discontinuities.

Yes, if the filtered signal still has discontinuities (but the width/duration of the "Gibbs ears" near the discontinuity points is zero). Discontinuities in the filtered signal can occur, for example, when a high-pass filter retains the edges of a square wave.

Of course, I'm referring to a continuous-time convolution and a continuous-time/continuous-frequency Fourier transform and inverse Fourier transform.

Your question does not apply to the discrete-time domain, since an ideal square wave is not bandwidth-limited and cannot be represented in the discrete-time domain. If you sample a square wave "as is", then you suffer from aliasing and a unique exact reconstruction from the samples is not possible. And if you apply a (continuous-time) anti-aliasing filter to the square wave before sampling, then you are no longer sampling a square wave, and the discrete-time FIR filter is not applied to a square wave either.

[gf]

So method (b) requires using appropriate windows.
Note that unless you need to filter only a relatively short buffer (and you don't need to do it in real time), using a single direct and inverse FT on the whole signal becomes unpractical and the usual way of doing it is the overlap-add method (and in particular the WOLA), in which you process a sequence of overlapping sections of signal:

https://www.dsprelated.com/freebooks/sasp/Overlap_Add_OLA_STFT_Processing.html

In order to do a simple FFT convolution with overlap-add ¹⁾, you must not use a synthesis window (WOLA), and a specific (non-rectangular) analysis window function is not required either. Only if time-dependent filtering or non-linear modifications in the frequency domain are desired, then an appropriate window function and/or WOLA are a good idea.

¹⁾ which is mathematically equivalent to a linear convolution of a (constant) FIR kernel with a (possibly infinite) stream of samples

[ By the way, with the STFT we are of course no longer working in the pure time domain or in the pure frequency domain, but in the time-frequency domain (i.e. we don't have a single spectrum for the whole infinite signal, but we are working with a sequence of time-dependent short-time spectra). ]

[Nominal Animal]

Just for fun, I created a CC0-1.0/Public Domain tool page, here, that works even if copied to your own computer, with no internet access), to experiment with the frequency spectrum needed to represent pulse trains, and the "Gibbsy"-like overshoot issues related to filters constructed using frequency response truncation.

The base pulse train is described by the infinite Fourier series with \$b_n = 4/(\pi\,n)\$ for odd \$n \ge 1\$, \$b_n = 0\$ for even \$n\$, and \$a_n = 0\$ for all \$n\$.  You can supply different constants in the input boxes, one per line, and when you tab away or click on the Calculate at top, it'll redraw the signal represented by the Fourier series coefficients.

At the background, you'll see a digital pulse with period 2.

If you use b_n = 4/(n*PI) 0 4/(n*PI) 0 4/(n*PI) 0 4/(n*PI) 0 4/(n*PI), a_n = 0 you can see what the background pulse would look like, when filtered by a low-pass filter that passes all frequencies at and below nine times the pulse frequency as-is, but blocks all frequencies at and above eleven times the pulse frequency completely.  The overshoot is the maximum, and pretty horrible.

The default 3.9/(n*PI) 0 3.4/(n*PI) 0 2.4/(n*PI) 0 1.5/(n*PI) 0 0.5/(n*PI) shows the effect of a theoretical filter that only affects amplitude, causes no phase shift, attenuating the pulse frequency by 1-3.9/4.0=2.5%, by 15% at triple the frequency, by 40% at five times the frequency, by 62.5% at seven times the frequency, by 87.5% at nine time the frequency, and blocks all frequencies above eleven times the pulse frequency.  The slew rate isn't high, but there is very little overshoot.

The idea here is to show that even if you had a perfect filter, to avoid overshoot you need to attenuate the higher frequencies at or above your signal frequency: you cannot have a steep stopband and no attenuation below it.  You could limit the slew rate instead –– corresponding to a filter whose frequency response varies depending on the signal.  And, that something like 10× the signal base frequency seems sufficient width for the attenuation band to recreate the original digital pulses (with very high input slew rate) with little to no overshoot, as an initial guess based on the default example.  I think you can construct a filter with narrower attenuation band by allowing phase distortion, but didn't delve into this that deep.
Finally, that devices with slew-rate limited outputs don't necessarily need more than say 10× the original digital signal frequency to be reproduced well, as shown by the frequency-limited default synthesized pulse shows.

[hamster_nz]

"The idea here is to show that even if you had a perfect filter, to avoid overshoot you need to attenuate the higher frequencies..."

Perhaps more concisely if you have a square wave (of any frequency) the range of fundamental will always overshoot the range of the square wave. Where the overshoots are, and by how much, depends what higher harmonics you filter, and the overall gain.

The 3rd harmonic pushes the overshoot to the towards the edges, the 5th harmonic does the same again. Eventually the overshoot ends up only at the edges.

You either have to have a gain less than 1 at your frequency of interest, or allow for overshoot.

(Graphic snipped from Wikipedia below)