[Fixed Thanks WP] Missing function in "ucomplex" "cabs"

[wp]

The line with cexp is wrong, the algorithm wants a complex argument but you are giving it a real argument.

After fixing this and plotting the original and filtered data in a chart, the noise reduction effect looks good (just changed the signal to something more interesting).

[marcov]

There are also pascal FFT routines for non power of two FFTs.  The site of Nils Haeck (also nativejpg and nativexml packages) seems gone, but there are still clones on github here and there.

[jamie]

Looking at the ucomplex overload for the ":=" and cexp and how the "cint(0,...) was used, one of them are backwards?

is it the ucomplex overload of ":=" or the C implementation?

Because the overload sets the Real, but you use the "cint" to set the imaginary number.

Jamie

[wp]

A complex number in uComplex is a record consisting of the real and imaginary parts. The normal assignment operator (inherited from the record type) assigns both real and imaginary parts. The overloaded special assignment operator of uComplex assigns only the real part, leaving the imaginary part zero:

Code: Pascal  [Select][+][-]

1.   operator := (r : real) z : complex;
2.     begin
3.        z.re:=r;
4.        z.im:=0.0;
5.     end;
6.  

What you need in the cexp call of the fft function, however, is a complex number with real part zero and non-zero imaginary part. The quickest way to get this is by using the cInit function:

Code: Pascal  [Select][+][-]

1. var
2.   z: complex;
3. begin
4.   z := cinit(0, -2.0*Pi*k/n);
5.   // or: z.re := 0.0;
6.   //     z.im := -2.0*pi*k/n;
7.   t := cExp(z) * ...;

[jamie]

Yes, I understand that, hence my question?

What is the gospel, ucomplex or the C implementation?

Jamie

[Nitorami]

In Pascal, complex numbers are most conveniently handled via a record. In C, "double complex" is a datatype from the C library <complex.h> and consists of a real and an imaginary part, both as double. I don't know how C does this internally, via a struct or not, but that does not matter. Both implementations are equivalent and do exactly the same.

Effectively what the C code does: -2.0 * I * M_PI * k / n is an angle i.e. a real.
cexp computes a complex re = cos(angle) and im = sin (angle).
In wp's example however, a complex is passed to cexp; cexp treats the imaginary part as angle, the real part as damping. That is an extension you don't need for the FFT, where damping must be zero.

[jamie]

I am having a problem conveying an issue here.

I understand what is happening under the hood.

Which one is wrong? Ucomplex or the C lib.

The math is correct eitherway, but the assignment to the fields is flipped flopped.

When I see an operator overload for ":=" to Memic an existing function, like in C, I would expect it to also assign the correct fields.

 Only one way works with this code, you can't have it either way.

 So, who has the correct implementation UComplex of ":=" overload or C?
Jamie

[Thaddy]

Yes, I understand that, hence my question?

What is the gospel, ucomplex or the C implementation?

Jamie

Both are the same:

Code: C  [Select][+][-]

1. double complex
2. /* this means a complex number with re and im as double */

Looks more like not understanding the C code.
The above is how C resolves the intended precision of the complex struct, which is an open typed macro in complex.h

[Nitorami]

Copilot says in C cexp (1) = e = 2.71 which means that C treats a real argument as real part i.e. as damping, same as in Pascal. But "I" is defined in C as the pure imaginary unit (re=0; im=1). In cexp(-2.0 * I * M_PI * k / n), multiplication with I would then deliver a complex number with re=0 and im as the angle -2*M_PI*k/n; same as in wp's example.

[Thaddy]

Posts crossed, but CoPilot is correct and therefor both Pascal and C. See my above reply. The error is misinterpreting the C code in the implementation.
wp already corrected that.

[jamie]

I have the C include file sitting Infront of me and I do understand C/C++ enough, and this is contradictory.

I'll conclude that no-one here really knows which one the correct action is to take on a ":="

This is getting out of hand.
but thanks to WP in any case for what he did.

[No more Please]

Jamie

[Thaddy]

Op. Cit.

Code: Text  [Select][+][-]

1. The line with cexp is wrong, the algorithm wants a complex argument but you are giving it a real argument.

No denials plz.

[jamie]

@Thaddy, I must say you are dense at times.

This has nothing to do with the cexp, it's all about the overloaded operator ":=" for the Complex record.

Its backwards which is why you can't directly enter the equation to the cexp, must build a Complex record outside that manually using cint or what ever.

do you get it now?

Jamie

[wp]

Found this: https://pubs.opengroup.org/onlinepubs/009695399/basedefs/complex.h.html

_Complex_I
    Expands to a constant expression of type const float _Complex, with the value of the imaginary unit (that is, a number i such that i2=-1).
imaginary
    Expands to _Imaginary.
_Imaginary_I
    Expands to a constant expression of type const float _Imaginary with the value of the imaginary unit.
I
    Expands to either _Imaginary_I or _Complex_I. If _Imaginary_I is not defined, I expands to _Complex_I.

In other words: The upper-case "I" in the expression "-2.0 * I * M_PI * k / n" is the "imaginary unit"; the term multiplied by I is the imaginary part. This way C has a very math-like way of defining complex variables: z = a + I * b (where a and b are floats). In FPC's uComplex, however, this possibility does not exist, we can only create complex numbers either by assigning the record elements one by one (z.re := a; z.im := b) or using the cInit function (z := cinit(a, b)).

If you'd declare a constant ImUnit: complex = (re:0.0, img:1.0) the expression could be written exactly like in C (it would not be very effective, though, since it introduces a lot of unnecessary multiplications and additions with 0):

Code: Pascal  [Select][+][-]

1. const
2.   ImUnit: complex = (re:0.0; im:1.0);
3. ...
4.   // Combine
5.   for k:= 0 to (n shr 1)-1 do
6.   begin
7.     t := cexp(-2.0 * ImUnit * pi * k / n)) * odd[k];  
8.     // I'd combine the real values, though: cexp(ImUnit*(-2.0 * pi * k / n)) * odd[k], but leaving this for clarity
9.     x[k] := even[k] + t;
10.     x[k + n div 2] := even[k] - t;
11.   end;
12.  

[EDIT]
No - I was wrong: in uComplex there is even a declaration for "i: complex = (re: 0.0; im:1.0)"; this way the C expression can be used literally. And there is a third way to define a complex number using uComplex:

Code: Pascal  [Select][+][-]

1. var
2.   z: complex;
3. ...
4.   z := 5.0 + 6.0*i;

[jamie]

And that fixes the problem.

I can use the code exactly the way it was intended now.

I use nameSpace uComplex.I  in the equation.

t := cexp(-2.0 * Ucomplex.I*PI * k / n) * odd[k];

To recap on the Complex := SomeFloatNumber, it is reversed to what the C Code uses but seems to work here.

Thank You very much. Maybe better Name for it would have been nice.