Cesaro summation of Fourier series to avoid Gibbs ripples

[Music Engineer]

Last night, I stumbled upon an interesting video from Trefor Bazett (the interesting part starts at 9:28 and ends at 11:28):

https://www.youtube.com/watch?v=AkPZcS8eqmA

and learned a potentially useful new bit of knowledge from it: How Cesaro summation applied to a Fourier series can be used to reduce the Gibbs ripples that occur when truncating the series - as we do when generating bandlimited waveforms. I knew about Cesaro summation before (from a Mathologer video) and, of course, I also knew about Fourier series before - but so far, I didn't make the connection. I immediately tried to implement the general idea - and, to my pleasant surprise, it seemed to work on the first try. The code is in my research codebase here:

https://github.com/RobinSchmidt/RS-MET- ... iments.cpp

in the function testFejerSum() (at line 9133 as of 2024/02/28). The code has been used to render these plots:

FejerSquare.png

FejerSquareZoomed.png

FejerSaw.png

I wonder, if this could be useful for audio processing - for example for generating (non-downsampled) mip-mapped wavetables via the usual FFT -> truncation -> IFFT approach. How would that Cesaro/Fejer approach compare to more standard techniques like tapering off the higher harmonics rather than doing a hard truncation? Did anyone ever do such a thing? Is this exciting or boring? Any other comments on this?

[Jeff McClintock]

i would be interested to see how the harmonic spectrum compares. i.e. does this technique reduce the higher partials?

[Music Engineer]

Yes - indeed. That's something I'd like to learn more about as well. I mean, the usual way to generate FFT/IFFT mip-maps is just spectral truncation (at least, that's the way, I'm currently doing it). One possible different way would be to apply some spectral tapering function - basically a sort of flat-top window on the spectrum before doing the IFFT...but this Cesaro averaging seems to be something very different [Edit: Nope - it doesn't - it applies a triangular window - see below]. I have not yet done much experimentation with that, so far - so I don't really know how it compares with respect to reduction of higher partials. But yeah - I would of course expect some sort of lowpass or high-shelf (attenuation) effect.

[juha_p]

Few years ago I did some comparisons between tapers in square wave generation. Here's one with couple tapers compared (Fejér and combined Lanczos-Fejér). Dunno if the plot there is correctly done :

Code: Select all

% needed for Octave ------------------------------    
pkg load signal
% ------------------------------------------------
clear all;
clf;

fs = 44100;
fc = 440;
w = 2 * pi * fc;
N = (fs/2)/fc;
r_start=0;
r_end=1;                            % length
r_step=1/fs;
t = r_start:r_step:r_end;           % range
r=t;
w_sqr = square(w*t);          % Square wave [-1:1]

figure(1);
plot(t,w_sqr);
axis([r_start r_end -1.2 1.2]);
grid on;
hold on;
audiowrite('testsquarewave.wav',w_sqr,fs);

% Approximation: Combined 
% Lanczos-Fejér tapers (Tukey window?) -----------
%n=(1:2:N);
%sumt_ = sum(((cos(pi.*n).-1).*(n.*2.*pi.*n.+(N.+1).*(-(2.*N.+1).*sin((2.*pi.*n)./...
%              (2.*N.+1)).-2.*pi.*n)).*sin(n.*w.*(0:1/fs:1)'))./((n.*(4.*pi^2.*N.+4.*pi^2).*n)),2);
              
% Fejér taper
n=(1:2:N)
sumt_ = sum((1-n/(N+1)).*(-1./(pi*n)).*sin(n.*w.*t').*(cos(pi*n)-1), 2);
              
Ffej = 2*sumt_;
Ffej=Ffej';

plot(t,Ffej);
axis([r_start r_end -1.2 1.2]);
grid on;
legend('Square', 'Fejér taper');
audiowrite('testFejer_squarewave.wav',Ffej,fs);

figure(2);
NFFT = 2*4096;
L=length(w_sqr);
X = fft(w_sqr,NFFT);
X = X(1:NFFT/2+1); %Throw the samples after NFFT/2 for single sided plot
Px1 = X.*conj(X)/(NFFT*L);
f = fs*(0:NFFT/2)/NFFT; %Frequency Vector
plot(f,10*log10(Px1),'g', 'linewidth', 1);
hold on;

L=length(Ffej);
X = fft(Ffej,NFFT);
X = X(1:NFFT/2+1); %Throw the samples after NFFT/2 for single sided plot
Px2 = X.*conj(X)/(NFFT*L);
f = fs*(0:NFFT/2)/NFFT; %Frequency Vector
plot(f,10*log10(Px2),'k','linestyle','--', 'linewidth', 1);

title('Single Sided Power Spectral Density');
xlabel('Frequency (Hz)')
ylabel('Power Spectral Density- P_{xx} dB/Hz');
legend('Square', 'Fejér taper');
ylim([-60 0]);

squarevsfejer.png

Noticed that one can "approximate" square wave and Fejér with an equation (add below code to above code):

Code: Select all

n=(1:2:N);
G = 2;  % G=2 => [-1,1] , G=1 => [-0.5,0.5]
a = 0.02;  % a=0.02 => close to Fejér, a=0.000001 => close to square wave
sumt_ = sum(G.*atan(sin(w.*t')/a)/pi + G/100000,2);
plot(t,sumt_, 'r', 'linestyle', '--');
axis([r_start r_end -1.2 1.2]);
grid on;
legend('Square', 'Fejér taper', '(G atan(sin(x)/a))/pi + G/10000');
audiowrite('atansin_squarewave.wav',sumt_,fs);

test02.png

test2.png

[signalsmith]

i would be interested to see how the harmonic spectrum compares. i.e. does this technique reduce the higher partials?

Yeah - that the Cesàro sum is an infinite limit of a sequence, but this is just showing the Nth value of that sequence, in which the i'th term in the sum is multiplied by max(0, 1 - i/N).

As N -> Infinity, the multiplier of each term -> 1, but if you're just taking the 15th element of the sequence, you get a linearly-ramped-down weighting of the first 15 harmonics. If you always matched N to Nyquist, the effective filter response looks like this.

[Music Engineer]

Ahhh...now I see! I had some misconception in my head about how this averaging of the different bandlimited time-domain signals translates to the effects on the spectrum. It actually translates to a spectral *windowing* or weighting - so doing this Cesaro summation thing is not actually doing anything new. It's just a different perspective on a known technique - namely, using spectral tapering windows to attenuate high frequency components. I will update my code soon with comments....

[Tone2 Synthesizers]

I do not think that Cesaro is a very good solution to get rid of the ripple. There seems to be too much loss in the high frequencies.
What I did to get rid of most ripple is to linear fade out the highest half/quarter of the harmonics.
What could be interesting to examine would be if/how phase shifts can further enhance the situation with the ripple

[Music Engineer]

Yeah - I now think that Cesaro summation just corresponds to using a triangular window (which is also known as the Fejer-window - probably that's why) in the frequency domain - which indeed implies a strong attenuation of high frequencies. Trying to use phase-shifting to counteract the ripple is indeed a thing that might be interesting to explore.

[Tone2 Synthesizers]

An idea could be to do something similar what a Bessel Filter does within the FFT transform

[Tone2 Synthesizers]

This looks good:
https://math.stackexchange.com/question ... phenomenon

[Tone2 Synthesizers]

I made an interesting observation: Unlike the Dirichlet the Fejer Kernel (that looks very similar to the coeficients of a FIR lowpass) does not have any negative values. I think this is the secret why it does not 'overshoot' and has only very few ringing.
https://en.wikipedia.org/wiki/Fej%C3%A9r_kernel

Does anyone know a closed form of the Fejer Kernel that can be directly applied to the fft frequency spectrum?

[Music Engineer]

Does anyone know a closed form of the Fejer Kernel that can be directly applied to the fft frequency spectrum?

If I get it right, the time domain "Dirichlet kernel" is a periodic version of the sinc function:

https://math.stackexchange.com/question ... let-kernel

The latter is the (inverse) Fourier transform of the rectangular window in the frequency domain. Likewise, it seems like the Fejer kernel is the periodic version of of the sinc^2 function - which corresponds to a triangular window function in the spectral domain.

Here, it says that the triangular window is also known as "Fejer window" (and also as "Bartlett window"):

https://en.wikipedia.org/wiki/Window_fu ... lar_window

...take that with a grain of salt, though.

However, I actually tried to synthesize my "Cesaro summed" output signal in the experiment above directly by applying a triangular window to the Fourier coeffs and it worked. The code has been updated. See the synthesis of "fejerWave2" - it is identical to "fejerWave" up to roundoff errors. See also the comment section "Conclusions" at the bottom of the function.

[mystran]

Does anyone know a closed form of the Fejer Kernel that can be directly applied to the fft frequency spectrum?

If I get it right, the time domain "Dirichlet kernel" is a periodic version of the sinc function:

https://math.stackexchange.com/question ... let-kernel

The latter is the (inverse) Fourier transform of the rectangular window in the frequency domain.

You obtain cardinal sine as you take the integral cos(pi*w*t)/2 dw over [-1,1]. Go ahead, work out the integral with pen and paper, it's not difficult (hint: use exp(i*w)=cos(w)+i*sin(w), saves the trouble of having to bother with trigonometry), it's quite elementary. This is literally the inverse Fourier transform of a rectangular function.

Similarly, you obtain a Dirichelet kernel when rather than taking an integral, you take a series... and that's the inserve Fourier series. The Fourier series is the periodic version of a Fourier transform in continuous domain the same as a discrete Fourier transform is a periodic version of a discrete-time (or more generally "domain") Fourier transform.

In other words.. we may take the time (or space) and spectrum to be either continuous or discrete independently of each other. Discrete spectra correspond to periodic waveforms. If you're feeling bored, you can also work out a proof that the Fourier transform of a periodic waveform must be zero except at isolated points corresponding to multiples of the base period.. and these isolated points are precisely the Fourier series. The cardinal sine and Dirichlet kernels are literally just (perhaps inverse, but in a sense all Fourier transforms are 90 degree rotations on the time-frequency plane so this is kinda irrelevant) transforms of a boxcar.

[mystran]

As for Fejér kernel, quoting https://en.wikipedia.org/wiki/Fej%C3%A9r_kernel "The traditional definition expresses the Fejér kernel F_n(x) in terms of the Dirichlet kernel..." with the expressions for nth Fejér kernel the average of the first n Dirichlet kernels. Now Dirichlet kernel is a boxcar (well, strictly an impulse train I guess) in spectral domain, so the average of progressively longer Dirichlet kernels has a triangular shape in terms of linear magnitude.

That's about as much as we need to care as far as audio goes. The fact that the kernel is non-negative everywhere could be of some use for some purposes (edit: I guess especially if you need a discrete summation formula), but it's periodic which limits it's applicability and in terms of actual audio signals the frequency response is frankly just garbage... and in terms of something like a square wave, it's not even particularly effective at reducing Gibb's other than pulling it down so that it doesn't overshoot... but it still oscillates anyway, so what's the point exactly?

Now, this is a bit opionated, but I think people obsess about Gibb's oscillations far more than is healthy... but if you actually must care, then why not just filter by a critically damped two-pole (which can be linear-phase if done in spectral domain for periodic waveforms)? Those have a strictly positive IR as well.. which doesn't oscillate the way Fejér kernel's do.

[Tone2 Synthesizers]

I agree with you that the frequency response is garbage. It is caused by the too-wide main lobe in the center of the kernel (and the missing negative values?). Making the main lobe more tight should fix the frequency response at least a bit.

Preventing overshoting makes a lot of sense if you process stuff with integers that can overflow.

My experience is that the gibbs phenomen can (but must not) have a negative effect on the sound quality.
It can make cubic splines used in following processing blocks sound worse. Splines do not like signals with lots of oscillation close to nyquist. I experienced a negative impact on the sound quality of some DA converters (which use splines for resampling?). I am pretty sure it can also have a negative impact on the sound of speakers. Speakers never are perfectly linear.
An oscillator with lots of gibbs ringing frequently has an audible impact on the sound of non-linear filters or anything else non-linear like distortions in the following processing blocks.
You also don't want to see such a ringing if you do accurate analog emulations. Analog oscillators simply do not have 'gibbs ringing'.

I do not think that any kind of ringing must be prevented at any cost. I think a healthy compromise is the way to go: Reducing the ringing without messing up the frequency response too much.