Fourier Transform in one BBC BASIC statement

[Hated Moron]

Would it surprise you to learn that you can implement a complete Fourier Transform calculation in one line of BBC BASIC code? Perhaps not (lines can, after all, be up to 251 bytes long). In that case, what if I was to say that you can implement a Fourier Transform in one BBC BASIC statement, would that impress you?

Well you can, and here is the code to prove it (in this case for an aperture of 1024 samples)! Firstly we need some initialisation, that's done just once so it doesn't count towards the number of statements per Fourier Transform (do sufficiently many and the contribution to the average from the initialisation tends to zero):

Code: Select all

      Aperture% = 1024
      DIM In(Aperture%-1), Out(Aperture%-1), DFT(Aperture%-1, Aperture%-1)
      REM. Prepare DFT coefficients:
      FOR I% = 0 TO Aperture%-2 STEP 2
        r = -PI*I%/Aperture%
        FOR J% = 0 TO Aperture%-1
          DFT(I%+0,J%) = COS(r*J%)
          DFT(I%+1,J%) = SIN(r*J%)
        NEXT
      NEXT I%

Now for the Fourier Transform calculation itself. Put 1024 input samples in the array In() then do this:

Code: Select all

      Out() = DFT() . In() 

The array Out() will contain the 512 real terms of the transform in the even-numbered elements and the 512 imaginary terms in the odd-numbered elements. Job done!

Well almost. In practice you will probably want to apply a window function (e.g. Hanning) to the input data, and if you want the magnitudes of the spectral components you will need to calculate SQR(real^2 + imaginary^2) but these are strictly not part of the Fourier Transform.

[Hated Moron]

On 05/08/2023 11:41, Kit Pearce via groups.io wrote (cross-posted from the Discussion Group):

Hi, any chance would you be able to provide code for the inverse transform indicating where we can zero frequencies to be removed?

So for the inverse-transform you need a complex input, yes? I did have that at one point but as a result of simplification I got rid of the (always zero) imaginary input terms, and sadly I didn't keep the intermediate version.

If you're at all mathematically adept I'm sure you can work it out from first principles. Once expressed as a sum of sines and cosines converting it to a dot-product is fairly straightforward.

One thing I would add is that if you're using it to process audio (in your case to notch out some frequency bands) the window function will introduce artefacts. To avoid those you would probably need to run two paths in parallel, one offset from the other by half the period of the window, and then sum again at the end.

Given that by the time you've done that, plus computing both forward and inverse transforms, you've increased the CPU demand by a factor of four, you probably wouldn't be able to handle real-time audio in BASIC unless you decrease the aperture and/or reduce the sampling rate significantly.

What exactly is it that you want to do? If it's simply to notch out a few frequencies, personally I'd be inclined to synthesise a FIR filter to do it; simpler and less computationally demanding than converting to the frequency domain and back. If you can tell me the specification you are trying to meet, I can attempt to synthesise a filter for you using FIRRBBC.

[Hated Moron]

personally I'd be inclined to synthesise a FIR filter to do it; simpler and less computationally demanding than converting to the frequency domain and back. If you can tell me the specification you are trying to meet, I can attempt to synthesise a filter for you using FIRBBC.

Just experimenting with FIR filter implementations in BBC BASIC (you can do that with a dot-product as well, of course) it certainly seems that a reasonable size filter is feasible at audio sampling frequencies. So it all depends on the performance required, most crucially how narrow a notch are you looking for? Do you need more than one?

[Hated Moron]

On 05/08/2023 11:41, Kit Pearce via groups.io wrote (cross-posted from the Discussion Group):

Hi, any chance would you be able to provide code for the inverse transform indicating where we can zero frequencies to be removed?

What exactly is it that you want to do? If it's simply to notch out a few frequencies, personally I'd be inclined to synthesise a FIR filter to do it; simpler and less computationally demanding than converting to the frequency domain and back. If you can tell me the specification you are trying to meet, I can attempt to synthesise a filter for you using FIRBBC.

As so often seems to happen, the OP (Kit Pearce) didn't respond to this. I'm still happy to help.