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

Hi,

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

Best regards,
Žilvinas

MSEgui has tfft which uses the FFTW library.
https://gitlab.com/mseide-msegui/mseide-msegui/blob/master/lib/common/math/msefft.pas

Martin

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

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

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

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...

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

-----
Many thanks ;-)

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.

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

-----
Many thanks ;-)

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.

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'

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

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

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.

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

-----
Many thanks ;-)

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

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 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