Good FFT example anywhere?

[Bo Berglund]

I am looking for some good example of FFT functions in pascal but it
seems like what Google serves up is rather old and refers to
Turbo-pascal and the like...
So maybe someone here knows about some open-source example of FFT
using FreePascal (or Delphi)?

I want to analyze the frequency content of transient responses
measured using a 24 bit A/D converter. It will produce 8192 samples
for each measurement.

Any suggestions welcome!

PS: I have saved this list from Sept 2003 and I searched it for the
word FFT without success, except the hits on the word offtopic.. DS


--
Bo Berglund
Developer in Sweden

_______________________________________________
fpc-pascal maillist  -  [hidden email]
http://lists.freepascal.org/cgi-bin/mailman/listinfo/fpc-pascal

[Žilvinas Ledas-2]

On 2017-04-09 11:06, Bo Berglund wrote:
> I am looking for some good example of FFT functions in pascal but it
> seems like what Google serves up is rather old and refers to
> Turbo-pascal and the like...
> So maybe someone here knows about some open-source example of FFT
> using FreePascal (or Delphi)?

Hi,

I was using attached unit some 7 years ago for various speech processing
purposes.


Best regards,
Žilvinas

>
> I want to analyze the frequency content of transient responses
> measured using a 24 bit A/D converter. It will produce 8192 samples
> for each measurement.
>
> Any suggestions welcome!
>
> PS: I have saved this list from Sept 2003 and I searched it for the
> word FFT without success, except the hits on the word offtopic.. DS
>
>

_______________________________________________
fpc-pascal maillist  -  [hidden email]
http://lists.freepascal.org/cgi-bin/mailman/listinfo/fpc-pascal

[Martin Schreiber-2]

On Sunday 09 April 2017 10:06:41 Bo Berglund wrote:
> I am looking for some good example of FFT functions in pascal but it
> seems like what Google serves up is rather old and refers to
> Turbo-pascal and the like...
> So maybe someone here knows about some open-source example of FFT
> using FreePascal (or Delphi)?
>
MSEgui has tfft which uses the FFTW library.
https://gitlab.com/mseide-msegui/mseide-msegui/blob/master/lib/common/math/msefft.pas

Martin
_______________________________________________
fpc-pascal maillist  -  [hidden email]
http://lists.freepascal.org/cgi-bin/mailman/listinfo/fpc-pascal

[Adriaan van Os-2]

Bo Berglund wrote:

> I am looking for some good example of FFT functions in pascal but it
> seems like what Google serves up is rather old and refers to
> Turbo-pascal and the like...
> So maybe someone here knows about some open-source example of FFT
> using FreePascal (or Delphi)?
>
> I want to analyze the frequency content of transient responses
> measured using a 24 bit A/D converter. It will produce 8192 samples
> for each measurement.
>
> Any suggestions welcome!

For maximum speed you would use vector-code. This could be assembly from Pascal. But it is easier
to use Intel IPP <https://software.intel.com/en-us/intel-ipp> or Apple Accelerate Framework
<https://developer.apple.com/reference/accelerate> when you are on OS X or iOS. Apple also has a
lib with C-source code <http://raquo.net/MatrixFFT-12.tar> that is much faster than the FFT in the
Accelerate Framework when the FFT size is large, see
<https://web.archive.org/web/20110806070107/http://images.apple.com/acg/pdf/FFTapps_20090909.pdf>.

Pascal interfaces for IPP are here <http://adriaan.biz/intel/ipp.pas.zip>.

There are also the FFTs in the GNU Scientific library
<https://www.gnu.org/software/gsl/manual/html_node/Fast-Fourier-Transforms.html> but it has license
restrictions <https://www.gnu.org/software/gsl/manual/html_node/GNU-General-Public-License.html>.

I have a Pascal FastFourierTransform.pas unit (without vector code), but it is not written by me.
You would have to contact the author for permission to use it. I can give the email address by
private email.

Regards,

Adriaan van Os

_______________________________________________
fpc-pascal maillist  -  [hidden email]
http://lists.freepascal.org/cgi-bin/mailman/listinfo/fpc-pascal

[Marco van de Voort]

In our previous episode, Bo Berglund said:
> Turbo-pascal and the like...
> So maybe someone here knows about some open-source example of FFT
> using FreePascal (or Delphi)?

When I looked around I was also pointed towards IPP and FFTW, but shied away
from the license consequences (reapllying every 12 months
https://software.intel.com/en-us/faq/free-software and GPL resp)

There is a fairly decent performing unit set under MPL out there by Nils
Haeck, that also performance well for samples sizes that are not a power of
two:

http://www.simdesign.nl/fft.html
 
_______________________________________________
fpc-pascal maillist  -  [hidden email]
http://lists.freepascal.org/cgi-bin/mailman/listinfo/fpc-pascal

[Adriaan van Os-2]

> I am looking for some good example of FFT functions in pascal but it
> seems like what Google serves up is rather old and refers to
> Turbo-pascal and the like...

But note that even a perfect FFT will suffer from some side-effects

<https://en.wikipedia.org/wiki/Gibbs_phenomenon>
<https://en.wikipedia.org/wiki/Ringing_artifacts>
<https://en.wikipedia.org/wiki/Window_function>

Regards,

Adriaan van Os
_______________________________________________
fpc-pascal maillist  -  [hidden email]
http://lists.freepascal.org/cgi-bin/mailman/listinfo/fpc-pascal

[Bo Berglund]

On Sun, 9 Apr 2017 11:22:00 +0300, Žilvinas Ledas
<[hidden email]> wrote:

>
>On 2017-04-09 11:06, Bo Berglund wrote:
>> I am looking for some good example of FFT functions in pascal but it
>> seems like what Google serves up is rather old and refers to
>> Turbo-pascal and the like...
>> So maybe someone here knows about some open-source example of FFT
>> using FreePascal (or Delphi)?
>
>I was using attached unit some 7 years ago for various speech processing
>purposes.
>
>Best regards,
>Žilvinas
>

I looked at it briefly and there are parts of it I don't understand,
for example:
NumberOfBitsNeeded, what purpose is this?
ReverseBits, does what and why?

I might have to run the code in a test program to see what is actually
going on. I did take the courses on Fourier transforms and the FFT
variant but that was a long time ago (like 40 years or so) so I am
rusty in that regard.

Generally (in reply to the other suggestions) I want a pure pascal
function so I can build it for the target platform, which will
probably end up being ARM and Linux (RPi3).
So no Intel stuff or Windows dll solutions...


--
Bo Berglund
Developer in Sweden

_______________________________________________
fpc-pascal maillist  -  [hidden email]
http://lists.freepascal.org/cgi-bin/mailman/listinfo/fpc-pascal

[fredvs]

Hello.

> I was using attached unit some 7 years ago for various speech processing purposes.

+ @Martin

Can those units be used for sound samples (array of float32 or array ot integer 32/16) ?

Thanks.

Fre;D

[Bo Berglund]

On Sun, 9 Apr 2017 08:01:07 -0700 (MST), fredvs
<[hidden email]> wrote:

>Hello.
>
>> I was using attached unit some 7 years ago for various speech processing
>> purposes.
>
>+ @Martin
>
>Can those units be used for sound samples (array of float32 or array ot
>integer 32/16) ?
>

The unit declares the fft procedure like this:

(*---------------------------------------------------------------------------
procedure fft

Calculates the Fast Fourier Transform of the array of complex numbers
represented by 'RealIn' and 'ImagIn' to produce the output complex
numbers in 'RealOut' and 'ImagOut'.
---------------------------------------------------------------------------*)
procedure fft (
    NumSamples:   word;   { must be a positive integer power of 2 }
    var  RealIn:   array of double;
    var  ImagIn:   array of double;
    var  RealOut:  array of double;
    var  ImagOut:  array of double );

As you can see it uses double, which means you will be OK to supply
both single and integer values, just load the array with them and they
will be converted.

--
Bo Berglund
Developer in Sweden

_______________________________________________
fpc-pascal maillist  -  [hidden email]
http://lists.freepascal.org/cgi-bin/mailman/listinfo/fpc-pascal

[fredvs]

> As you can see it uses double, which means you will be OK to supply
> both single and integer values, just load the array with them and they will be converted.

Many thanks for answer.

Hmmm, it seems magic and very simple to use.
Ok, I will try with my audio library.

By the way (who knows) did you try it, for example, to isolate some frequency from a buffer of sample (like a high-pass filter) ?
(Of course example are welcome ;-) ).

Fre;D

[Rolf Grunsky]

On 2017-04-09 04:06 AM, Bo Berglund wrote:

> I am looking for some good example of FFT functions in pascal but it
> seems like what Google serves up is rather old and refers to
> Turbo-pascal and the like...
> So maybe someone here knows about some open-source example of FFT
> using FreePascal (or Delphi)?
>
> I want to analyze the frequency content of transient responses
> measured using a 24 bit A/D converter. It will produce 8192 samples
> for each measurement.
>
> Any suggestions welcome!
>
> PS: I have saved this list from Sept 2003 and I searched it for the
> word FFT without success, except the hits on the word offtopic.. DS
>
>

The book "Numerical Recipes: The Art of Scientific Computing" has an
entire section of FFT. The first edition of the book has code in both
FORTRAN and Pascal. You can find out more at "numerical.recipes". Later
editions are C and C++ and may not have Pascal examples.

There are many other scientific computing topics in the book.

--
                                TRUTH in her dress finds facts too tight.
                                In fiction she moves with ease.
                                Stray Birds by Rabindranath Tagore
_______________________________________________
fpc-pascal maillist  -  [hidden email]
http://lists.freepascal.org/cgi-bin/mailman/listinfo/fpc-pascal

[Ched]

Hi,

As an owner of the Numerical Receipes in Pascal, I can say that the sources codes are *very very old*
styled and a rewritting to dynamical array of data is a good idea ! The N.R. are probably written by
ForTranists :) But ths book contains lots of useful knowledge en tricks for numericians and is pleasant
to read.

Cheers, Ched'





Le 09.04.2017 à 19:21, Rolf Grunsky a écrit :

> On 2017-04-09 04:06 AM, Bo Berglund wrote:
>> I am looking for some good example of FFT functions in pascal but it
>> seems like what Google serves up is rather old and refers to
>> Turbo-pascal and the like...
>> So maybe someone here knows about some open-source example of FFT
>> using FreePascal (or Delphi)?
>>
>> I want to analyze the frequency content of transient responses
>> measured using a 24 bit A/D converter. It will produce 8192 samples
>> for each measurement.
>>
>> Any suggestions welcome!
>>
>> PS: I have saved this list from Sept 2003 and I searched it for the
>> word FFT without success, except the hits on the word offtopic.. DS
>>
>>
> The book "Numerical Recipes: The Art of Scientific Computing" has an entire section of FFT. The first
> edition of the book has code in both FORTRAN and Pascal. You can find out more at "numerical.recipes".
> Later editions are C and C++ and may not have Pascal examples.
>
> There are many other scientific computing topics in the book.
>

_______________________________________________
fpc-pascal maillist  -  [hidden email]
http://lists.freepascal.org/cgi-bin/mailman/listinfo/fpc-pascal

[Giuliano Colla]

Il 09/04/2017 10:06, Bo Berglund ha scritto:

> I am looking for some good example of FFT functions in pascal but it
> seems like what Google serves up is rather old and refers to
> Turbo-pascal and the like...
> So maybe someone here knows about some open-source example of FFT
> using FreePascal (or Delphi)?
>
> I want to analyze the frequency content of transient responses
> measured using a 24 bit A/D converter. It will produce 8192 samples
> for each measurement.
>
> Any suggestions welcome!
>
> PS: I have saved this list from Sept 2003 and I searched it for the
> word FFT without success, except the hits on the word offtopic.. DS
>
>

For mechanical vibration analysis (with a similar number of samples but
only 16 bit resolution) I have used the FFT library of Nils Haeck which
you may find here: http://www.simdesign.nl/fft.html
For my needs it was appropriate. You may give a look. It comes with
Mozilla Licence.

Regards,

Giuliano

_______________________________________________
fpc-pascal maillist  -  [hidden email]
http://lists.freepascal.org/cgi-bin/mailman/listinfo/fpc-pascal

[Bo Berglund]

On Sun, 9 Apr 2017 11:34:53 -0700 (MST), fredvs
<[hidden email]> wrote:

>> As you can see it uses double, which means you will be OK to supply
>> both single and integer values, just load the array with them and they
>> will be converted.
>
>Many thanks for answer.
>
>Hmmm, it seems magic and very simple to use.
>Ok, I will try with my audio library.
>
>By the way (who knows) did you try it, for example, to isolate some
>frequency from a buffer of sample (like a high-pass filter) ?
>(Of course example are welcome ;-) ).

Forgot to say that the code limits the number of samples like this:

procedure fft (
    NumSamples:   word;   { must be a positive integer power of 2 }

and

function IsPowerOfTwo ( x: word ): boolean;
var   i, y:  word;
begin
    y := 2;
    for i := 1 to 15 do begin <== Does not go to 16

So you are limited to 32768 since the next multiple of 2 (65536) is
too big to fit inside a word size variable.


--
Bo Berglund
Developer in Sweden

_______________________________________________
fpc-pascal maillist  -  [hidden email]
http://lists.freepascal.org/cgi-bin/mailman/listinfo/fpc-pascal

[fredvs]

> So you are limited to 32768 since the next multiple of 2 (65536) is too big to fit inside a word size variable.

Ha, ok, usually, by default, I use 65536 for all DSP's but even with 1024 I get acceptable result (and better latency).

OK, I will try.

Write you asap.

Fre;D

[Michael Schnell]

On 09.04.2017 10:06, Bo Berglund wrote:
> It will produce 8192 samples
> for each measurement.
If the "measurements" come as a sequence of sample-blocks, you
additional to the core FFT will need a sliding Window algorithm to e.g.
create something like a life spectrum display or do a useful convolution.

-Michael
_______________________________________________
fpc-pascal maillist  -  [hidden email]
http://lists.freepascal.org/cgi-bin/mailman/listinfo/fpc-pascal

[Bo Berglund]

On Mon, 10 Apr 2017 09:58 +0200, Michael Schnell
<[hidden email]> wrote:

>>It will produce 8192 samples
>>for each measurement.
>
> If the "measurements" come as a sequence of sample-blocks, you
> additional to the core FFT will need a sliding Window algorithm to e.g.
> create something like a life spectrum display or do a useful convolution.

Well, it is actually like this:
- A pulse stimulus is applied to the object to measure
- The response transient is recorded for up to 150 ms

The task now is to analyze the single transient for its constituent
frequencies (amplitude and phase) because these hold information about
the measured object.

The used 24 bit A/D converter can sample up to 50 KSMPS, but we don't
know if this is really necessary. In any case 164 ms of 50 kHz
sampling gives us 8192 samples to work with, reduce the rate to half
and we have 4096 samples for the same time instead.


--
Bo Berglund
Developer in Sweden

_______________________________________________
fpc-pascal maillist  -  [hidden email]
http://lists.freepascal.org/cgi-bin/mailman/listinfo/fpc-pascal

[Michael Schnell]

On 10.04.2017 23:20, Bo Berglund wrote:
> - A pulse stimulus is applied to the object to measure
> - The response transient is recorded for up to 150 ms

OK. In this case you will not need the "Windowing".

Instead you need to make sure that the pulse is surrounded by enough
Zero-line data. I suggest to place the Pulse's non-zero data in the
middle of an array of at least the double count of samples, the rest
being zero. Supposedly th Pulse starts with Zero, anyway, so no problem
here. If at the end the measurement is truncated, you should use a
"sinc" or similar curve to force it smoothly down to zero with the last
measured sample before the Zero-line starts.

Not doing this will result in a huge aliasing and unusable "spectrum data.

The cause is that Digital Fourier Transform converts a number of samples
in the same number of "bins" (= complex amplitude / phase spectrum data
for a dedicated frequency).

If there are n (="Window-Size") samples and the Sample Frequency is f,
the bins represent the frequencies 0, f/n, 2f/n, 3f/n, ... Hence the
bin[1] represents the Sample Frequency by the Window Size, which is the
lowest frequency usable, while the last bin represents the Sample
Frequency which in fact is not usable due to the Nyquist Theorem, that
allows for using only frequencies lower than half the Sample Frequency.
(This said, you need to make sure by hardware means that no frequencies
higher than half the sample even reach the A/D converter. Otherwise you
will feature aliasing that will make the spectrum unusable. )

The Nyquist Theorem also is visible in the fact that a DFT does not
really convert the signal you measured to a spectrum you want to see,
but in fact converts a periodic signal to a periodic spectrum. The
(assumed) period of the signal being the Window Size divided by the
Sample Frequency, the period of the spectrum being the Sample Frequency
divided by the Window Size. Hence adding the Zero-Line will allow your
signal to be the significant part of a period.

-Michael
_______________________________________________
fpc-pascal maillist  -  [hidden email]
http://lists.freepascal.org/cgi-bin/mailman/listinfo/fpc-pascal