Detect Frequency with FFT

Share your advanced PureBasic knowledge/code with the community.
chris319
Enthusiast
Enthusiast
Posts: 782
Joined: Mon Oct 24, 2005 1:05 pm

Detect Frequency with FFT

Post by chris319 »

Here is a small program which detects the dominant frequency in a wav file. The wav file must be a 16-bit monaural file. You can fool around with a sine wave for starters to see how it works.

Code: Select all

;OPENS A 16-BIT MONAURAL WAV FILE CONTAINING A SINE WAVE
;DETECTS THE FREQUENCY OF THE SINE WAVE
;READS SAMPLE RATE OF FILE FROM format\nSamplesPerSec
;RUNS IN CONSOLE MODE
;CREATED 8/9/08 by chris319


Global format.WAVEFORMATEX
Global *soundbuffer, bytes_read, real.d, imag.d, frad.d, frequency.d, imd.d, amplitude.d, max_amp.d, peak_freq.d
Global filename$
Global Dim SampleData.d(1)
#PI2 = 2 * #PI

OpenConsole()

Procedure load_wav()
filesize = FileSize(filename$) - 66
*soundbuffer = AllocateMemory(filesize)
ReadFile(1, filename$)
FileSeek(1, 20) ;SKIP HEADERS
ReadData(1, @format, SizeOf(format))
FileSeek(1, 68) ;SKIP HEADERS
bytes_read = ReadData(1, *soundbuffer, filesize)
CloseFile(1)
ReDim SampleData.d(bytes_read / 2)

ct2 = 0
For ct = 0 To bytes_read / 2

  SampleData(ct) = PeekW(*soundbuffer + ct2) / 32767 ;READ INTO ARRAY AND NORMALIZE SAMPLE T0 32767

  ct2 = ct2 + 2
Next

FreeMemory(*soundbuffer)
EndProcedure

Procedure FFT(frequency)
real.d = 0
imag.d = 0
frad.d = (frequency * #PI2) / format\nSamplesPerSec
samp_ct = 0

;'real' = real part of answer 
;'imag' = imaginary part of answer 
;'frequency' specifies the frequency you're looking for in radians/sec.

;HERE WE SWEEP THE SAMPLES TO FIGURE OUT THE AMPLITUDE AT THE SPECIFIED FREQUENCY
  For sample = 0 To format\nSamplesPerSec
    real = real + SampleData(sample) * Sin(frad * samp_ct)
    imag = imag + Sin(frad * samp_ct)
    samp_ct + 1
  Next

;Amount of 'frequency' present is then SquareRoot(real^2 + imag^2)
amplitude.d = Sqr(Pow(real, 2) + Pow(imag, 2))
EndProcedure

Print("Enter name of wav file: "): filename$ = Input(): load_wav()

CreateFile(1, "fft_out.csv")
WriteString(1, "Frequency,Amplitude" + Chr(10))

max_amp = 0

frequency = 100
While frequency <= 20000
  FFT(frequency)
  If amplitude > max_amp
    max_amp = amplitude
    peak_freq = frequency
  EndIf
  WriteString(1, StrD(frequency,0) + "," + StrD(amplitude,4) + Chr(10))
  frequency + 100
Wend

CloseFile(1)

PrintN (Chr(10) + "Peak frequency: " + StrD(peak_freq,0) + " Hz" + Chr(10) + Chr(10) + "Done")

here: Delay(250): Goto here

End
untune
User
User
Posts: 56
Joined: Sun Jul 06, 2008 10:07 pm
Location: Bolton, UK

Post by untune »

Hi Chris

This looks cool, I'm going to take a look at the tomorrow. I'm interested to see what PB can do with sound, when I finish my GUI/graphics system I'm hoping to tie them together in some way :)
KarLKoX
Enthusiast
Enthusiast
Posts: 681
Joined: Mon Oct 06, 2003 7:13 pm
Location: France
Contact:

Post by KarLKoX »

Really nice, i played a little bit with FFT too, i ve done something a long time ago using the F64 lib (when PB lacks 64 bits type, thanks to Jack ;) ), if you can read dirty code, read my code :P (FFT_F64() is the working func)
"Qui baise trop bouffe un poil." P. Desproges

http://karlkox.blogspot.com/
Froggerprogger
Enthusiast
Enthusiast
Posts: 423
Joined: Fri Apr 25, 2003 5:22 pm
Contact:

Post by Froggerprogger »

I cannot see the FFT here. You 'sweep' for 200 frequencies each time the whole sampledata for creating an amplitude-value for each. But the FFT defines an algorithm to get these values much faster by a concrete divide-and-conquer-algorithm that is not implemented here.
But anyway, this code is cool and works for arbitrary frequencies and sampledatalengths!
%1>>1+1*1/1-1!1|1&1<<$1=1
chris319
Enthusiast
Enthusiast
Posts: 782
Joined: Mon Oct 24, 2005 1:05 pm

Post by chris319 »

I cannot see the FFT here. You 'sweep' for 200 frequencies each time the whole sampledata for creating an amplitude-value for each. But the FFT defines an algorithm to get these values much faster by a concrete divide-and-conquer-algorithm that is not implemented here.
You are absolutely right. This is not a true FFT but I don't know what to call it. Is it a DFT (discrete Fourier transform)? Anyway, I am finding it very useful for audio measurements.
Froggerprogger
Enthusiast
Enthusiast
Posts: 423
Joined: Fri Apr 25, 2003 5:22 pm
Contact:

Post by Froggerprogger »

It is not even a FT :) The DFT gives the maths for presenting a signal as a sum of many discrete sinuses and cosinuses at some complex roots of unity. The FFT describes an algorithm to calculate the coefficients of the DFT in a fast way by divide-and-conquer that depends on some characteristics of these roots. But the frequencies under watch are fixed in the DFT/FFT and depend on the number of samples for the transformation, and cannot be chosen arbitrary.
%1>>1+1*1/1-1!1|1&1<<$1=1
Post Reply