Analyzing audio from specific source in real-time

[Boleeman]

Also found an interesting Delphi SinSyn - Sinusoidal Synthesizer here:

https://github.com/tothpaul/Delphi-SinSyn

Has HCF (Highest common factor) measurement

[Thaddy]

Good luck finding a pure "native" pascal solution.

Well, FredericVanMol's - Fruity Studio fame - FFT/IFFT is fast pure pascal and buffer scaling to display is very easy anyway.
I use it since the early 2000's for VST plugins and wav/mp3 players. And that is by no means the only pure pascal implementation. There are many more. One complete example is in the AsioVST package. Also the DSP* components contain a similar component.

[Thaddy]

Not a replacement for FFTW, but a simple one. Fast enough for display, even fast enough for audio processing if the buffers aren't too big.
Try to make use of the compiler optimizations. FPU choices help a lot.

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

1. unit tdkfft;
2. {
3.   Simple Cooley-Tukey (Gaussian) FFT/IFFT + scaling
4.   Fast enough for display buffers (256/512/1024/2048
5.   Also fast enough for audio processing if the buffers aren't too big.
6.   Simpler than the other options but very lightweight.
7.   All arrays must be N = a power of two
8.  
9.   Have fun,
10.  
11.   Thaddy
12. }{$mode objfpc}
13. interface
14. type
15.   // adjust here
16.   fft_t = double;
17.   TComplex = record
18.     Re, Im: fft_t;
19.   end;
20.  
21.   TComplexArray = array of TComplex;
22.  
23. // helpers
24. function Complex(const Re, Im: fft_t): TComplex;inline;
25. operator + (const a, b: TComplex): TComplex;inline;
26. operator - (const a, b: TComplex): TComplex;inline;
27. operator * (const a, b: TComplex): TComplex;inline;
28. function Scale(const a: TComplex; s: fft_t): TComplex;inline;
29. //FFT IFFT
30. procedure FFT(var Data: TComplexArray; N: Integer; Inverse: Boolean);inline;
31.  
32. implementation
33.  
34. function Complex(const Re, Im: fft_t): TComplex;inline;
35. begin
36.   Result.Re := Re;
37.   Result.Im := Im;
38. end;
39.  
40. operator + (const a, b: TComplex): TComplex;inline;
41. begin
42.   Result.Re := a.Re + b.Re;
43.   Result.Im := a.Im + b.Im;
44. end;
45.  
46. operator - (const a, b: TComplex): TComplex;inline;
47. begin
48.   Result.Re := a.Re - b.Re;
49.   Result.Im := a.Im - b.Im;
50. end;
51.  
52. operator * (const a, b: TComplex): TComplex;inline;
53. begin
54.   Result.Re := a.Re * b.Re - a.Im * b.Im;
55.   Result.Im := a.Re * b.Im + a.Im * b.Re;
56. end;
57.  
58. function Scale(const a: TComplex; s: fft_t): TComplex;inline;
59. begin
60.   Result.Re := a.Re * s;
61.   Result.Im := a.Im * s;
62. end;
63.  
64. procedure FFT(var Data: TComplexArray; N: Integer; Inverse: Boolean);inline;
65. var
66.   i, j, k, m: Integer;
67.   theta: Double;
68.   w, wp, temp: TComplex;
69. begin
70.   if N <= 1 then Exit;
71.  
72.   // Bit-reversal permutation
73.   j := 0;
74.   for i := 0 to N - 1 do
75.   begin
76.     if j > i then
77.     begin
78.       temp := Data[i];
79.       Data[i] := Data[j];
80.       Data[j] := temp;
81.     end;
82.     m := N div 2;
83.     while (m >= 1) and (j >= m) do
84.     begin
85.       j := j - m;
86.       m := m div 2;
87.     end;
88.     j := j + m;
89.   end;
90.  
91.   // Cooley-Tukey FFT
92.   m := 1;
93.   while m < N do
94.   begin
95.     theta := Pi / m;
96.     if Inverse then theta := -theta;
97.     wp := Complex(Cos(theta), Sin(theta)); // Correct twiddle factor
98.     w := Complex(1.0, 0.0);
99.     for k := 0 to m - 1 do
100.     begin
101.       i := k;
102.       while i < N do
103.       begin
104.         j := i + m;
105.         temp := w * Data[j];
106.         Data[j] := Data[i] - temp;
107.         Data[i] := Data[i] + temp;
108.         i := i + 2 * m;
109.       end;
110.       w := w * wp;
111.     end;
112.     m := m * 2;
113.   end;
114.  
115.   // Scale the result for IFFT
116.   if Inverse then
117.     for i := 0 to N - 1 do
118.       Data[i] := Scale(Data[i], 1.0 / N);
119. end;
120.  
121. end.

Demo:

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

1. {$mode objfpc}{$I-}
2. uses tdkfft;
3. procedure PrintComplexArray(const Data: TComplexArray; N: Integer);inline;
4. var
5.   i: Integer;
6. begin
7.   for i := 0 to N - 1 do
8.     WriteLn(Data[i].Re:1:10,' ', Data[i].Im:1:10);
9. end;
10.  
11. var
12.   Data: TComplexArray;
13.   N, i: Integer;
14. begin
15.   N := 8; // Must be a power of 2
16.   SetLength(Data, N);
17.  
18.   // Example input data (a simple signal)
19.   for i := 0 to N - 1 do
20.     Data[i] := Complex(Sin(2 * Pi * i / N), 0.0);
21.  
22.   WriteLn('Original Data:');
23.   PrintComplexArray(Data, N);
24.  
25.   // Compute FFT
26.   FFT(Data, N, False);
27.   WriteLn('FFT Result:');
28.   PrintComplexArray(Data, N);
29.  
30.   // Compute IFFT
31.   FFT(Data, N, True);
32.   WriteLn('IFFT Result (should match original data):');
33.   PrintComplexArray(Data, N);
34. end.