| Still working on the code for window functions. I've now discovered a bug in
the code for the sinc window. This is current code:
case 9: /* Sinc */
arg = TWOPI / flen;
for (i = -PI ; i < 0 ; i += arg)
*ft++ = (MYFLT)(xarg * sin(i) / i);
*ft++ = FL(1.0);
for (i = arg ; i < PI ; i += arg)
*ft++ = (MYFLT)(xarg * sin(i) / i);
return;
The reason why the *ft++ = FL(1.0); is used in between the two loops is that
the calculation ( sin(i) / i ) is illegal for i=0. From math calculus (using
L'Hopital's rule) it is known that
lim ( sin(i) / i) = 0 when i -> 0
In gen20 amplitude xarg is implemented, and hence this value has to be
scaled:
*ft++ = FL(1.0 * xarg);
The current bug can be seen clearly in Kevin's illustration of the various
window functions:
http://www.kevindumpscore.com/docs/csound-manual/miscwindows.html
In the sinc example, max amplitude = 0.75, but at i=0 there's an erraneous
spike hitting the value 1.
Corrected code should read:
case 9: /* Sinc */
arg = TWOPI / flen;
for (i = -PI ; i < 0 ; i += arg)
*ft++ = (MYFLT)(xarg * sin(i) / i);
*ft++ = FL(xarg); // <------- This line is changed
for (i = arg ; i < PI ; i += arg)
*ft++ = (MYFLT)(xarg * sin(i) / i);
return;
Best,
Trond
On 2004-08-18 10:41, "Trond Lossius" wrote:
> Hi!
>
> I'm currently porting the gen20 window functions from fgens.c to Java in order
> to use them in Max 4.5 inside the mxj object. While doing this I'm wondering
> about a detail of the implementations:
>
> In spite of the mathematical functions for the various windows being
> symetrical, the implementation is not. E.g. Hamming, Hanning, etc. are
> calculated this way:
>
> for (i = 0.0 ; i < TWOPI ; i += arg)
> *ft++ = (MYFLT)(xarg * (cf[0] - cf[1]*cos(i) +
> cf[2]*cos(2.0*i) - cf[3]*cos(3.0*i)));
>
> In this loop i will stop just short of TWOPI, and thus the last sample do not
> equal the first.
>
> The same goes for the Gaussian window:
>
> case 6: /* Gaussian */
> arg = 12.0/flen;
> for (i = -6.0 ; i < 0.0 ; i += arg)
> *ft++ = (MYFLT)(xarg * (pow( 2.71828, -(i*i)/2.0)));
> for (i = arg ; i < 6.0 ; i += arg)
> *ft++ = (MYFLT)(xarg * (pow( 2.71828, -(i*i)/2.0)));
> return;
>
> The counter will be running from -6 (inclusive) to +6 (exclusive)
>
> It probably doesn't make much of a difference for large buffers, but still...
>
> Best,
> Trond |