Get frequency in sound using _Vince algorithm

[Petr]

In the first place I want to thank @_vince  again for his program. This is just a modification of its code. I made from his example a function to make it as easy to use as possible. Use any stereo audio file supported in QB64 to try it out.

Code: QB64: [Select]

1.  
2. '''
3. DIM Left AS _MEM
4. DIM Right AS _MEM
5. DIM Analyzer(8) AS SINGLE  '   for left audio channel
6. DIM Analyzer2(8) AS SINGLE '  for right audio channel
7. '''
8.  
9. '''
10. AudioFile$ = "01.mp3"
11. S = _SNDOPEN(AudioFile$)
12. Left = _MEMSOUND(S, 1)
13. Right = _MEMSOUND(S, 2)
14.  
15. '''
16. SCREEN _NEWIMAGE(1024, 800, 32)
17. _SNDPLAY S
18. DO UNTIL _SNDPLAYING(S) = 0
19.     CLS
20.     FL = GetFreQ(S, Left, 1024)
21.     FR = GetFreQ(S, Right, 1024)
22.     _PRINTSTRING (0, 50), "Frequention on left:" + STR$(FL)
23.     _PRINTSTRING (0, 70), "Frequention on right:" + STR$(FR)
24.  
25.  
26.     SELECT CASE FL '                                                                              left channel
27.         CASE 20 TO 125: Analyzer(1) = Analyzer(1) + 30 'areas set by picture on google...
28.         CASE 126 TO 300: Analyzer(2) = Analyzer(2) + 30
29.         CASE 300 TO 500: Analyzer(3) = Analyzer(3) + 30
30.         CASE 500 TO 1000: Analyzer(4) = Analyzer(4) + 30
31.         CASE 1001 TO 2000: Analyzer(5) = Analyzer(5) + 30
32.         CASE 2001 TO 4000: Analyzer(6) = Analyzer(6) + 30
33.         CASE 4001 TO 8000: Analyzer(7) = Analyzer(7) + 30
34.         CASE IS > 8001: Analyzer(8) = Analyzer(8) + 30
35.     END SELECT
36.  
37.     SELECT CASE FR '                                                                              right channel
38.         CASE 20 TO 125: Analyzer2(1) = Analyzer2(1) + 30 'areas set by picture on google...
39.         CASE 126 TO 300: Analyzer2(2) = Analyzer2(2) + 30
40.         CASE 300 TO 500: Analyzer2(3) = Analyzer2(3) + 30
41.         CASE 500 TO 1000: Analyzer2(4) = Analyzer2(4) + 30
42.         CASE 1001 TO 2000: Analyzer2(5) = Analyzer2(5) + 30
43.         CASE 2001 TO 4000: Analyzer2(6) = Analyzer2(6) + 30
44.         CASE 4001 TO 8000: Analyzer2(7) = Analyzer2(7) + 30
45.         CASE IS > 8001: Analyzer2(8) = Analyzer2(8) + 30
46.     END SELECT
47.  
48.  
49.  
50.     'analyzer overflow control
51.     FOR Analyze = 1 TO 8 '                                    analyze 8 spectrums
52.         IF Analyzer(Analyze) > 150 THEN Analyzer(Analyze) = 150
53.         IF Analyzer2(Analyze) > 150 THEN Analyzer2(Analyze) = 150
54.     NEXT
55.  
56.     'draw analyzer
57.     AnC = 0
58.     FOR AnalyzatorX = 600 TO 670 STEP 10
59.         AnC = AnC + 1
60.         FreqHeight = Analyzer(AnC) '                          Y coordinate for drawing frequency analyzer, maximal size is 150 pixels
61.         FreqHeight2 = Analyzer2(AnC)
62.         LINE (AnalyzatorX - 4, 400)-(AnalyzatorX + 4, 400 - FreqHeight), &HFFAD0000, BF
63.         LINE (AnalyzatorX + 300 - 4, 400)-(AnalyzatorX + 300 + 4, 400 - FreqHeight2), &HFFAD0000, BF
64.         IF Analyzer(AnC) > 0 THEN Analyzer(AnC) = Analyzer(AnC) - 1 '
65.         IF Analyzer2(AnC) > 0 THEN Analyzer2(AnC) = Analyzer2(AnC) - 1 ' always, if analyzer contains some higher position than bottom, give one pixel out
66.     NEXT
67.  
68.     _DISPLAY
69.     _LIMIT 200
70. LOOP
71.  
72. SYSTEM
73.  
74. FUNCTION GetFreQ (SoundHandle AS LONG, arr AS _MEM, Samples AS LONG)
75.     sw = Samples&
76.     DIM c AS DOUBLE
77.     DIM u_r AS DOUBLE, u_i AS DOUBLE
78.     DIM v_r AS DOUBLE, v_i AS DOUBLE
79.  
80.  
81.     DIM pi AS DOUBLE
82.     pi = 4 * ATN(1)
83.     DIM x_r(sw - 1), x_i(sw - 1)
84.     DIM xx_r(sw - 1), xx_i(sw - 1)
85.     DIM t AS DOUBLE
86.  
87.     Position& = _SNDGETPOS(SoundHandle) * _SNDRATE * 2 '                 *2, because sound data are in INTEGER format and INTEGER is always 2 byte long
88.     IF Position& > Samples& AND Position& < arr.SIZE - Samples& THEN 'current (actual playing) sound frame is in middle - in position Position&
89.         j& = Position& - Samples&
90.         DO UNTIL i >= Samples&
91.             x_r(i) = 100 * _MEMGET(arr, arr.OFFSET + j&, INTEGER) / 32767
92.             i = i + 1
93.             j& = j& + 2
94.         LOOP
95.         '        test& = j& - Position&
96.     END IF
97.     i = 0
98.     j& = 0
99.     rfft xx_r(), xx_i(), x_r(), sw
100.     'find peak
101.     DIM max AS DOUBLE, d AS DOUBLE
102.     max = 0
103.     m = 0
104.     FOR i = 0 TO sw / 2
105.         d = 0.01 * SQR(xx_r(i) * xx_r(i) + xx_i(i) * xx_i(i))
106.         IF d > max THEN
107.             max = d
108.             m = i
109.         END IF
110.     NEXT
111.     IF m = 0 THEN GOTO sem
112.     u_r = xx_r(m - 1) - xx_r(m + 1)
113.     u_i = xx_i(m - 1) - xx_i(m + 1)
114.     v_r = 2 * xx_r(m) - xx_r(m - 1) - xx_r(m + 1)
115.     v_i = 2 * xx_i(m) - xx_i(m - 1) - xx_i(m + 1)
116.     c = (u_r * v_r + u_i * v_i) / (v_r * v_r + v_i * v_i)
117.     sem:
118.     GetFreQ = (m + c) * _SNDRATE / Samples& '                                                             1024 = number of samples
119. END FUNCTION
120.  
121. SUB rfft (xx_r(), xx_i(), x_r(), n)
122.     DIM w_r AS DOUBLE, w_i AS DOUBLE, wm_r AS DOUBLE, wm_i AS DOUBLE
123.     DIM u_r AS DOUBLE, u_i AS DOUBLE, v_r AS DOUBLE, v_i AS DOUBLE
124.  
125.     log2n = LOG(n / 2) / LOG(2)
126.  
127.     FOR i = 0 TO n / 2 - 1
128.         rev = 0
129.         FOR j = 0 TO log2n - 1
130.             IF i AND (2 ^ j) THEN rev = rev + (2 ^ (log2n - 1 - j))
131.         NEXT
132.  
133.         xx_r(i) = x_r(2 * rev)
134.         xx_i(i) = x_r(2 * rev + 1)
135.     NEXT
136.  
137.     FOR i = 1 TO log2n
138.         m = 2 ^ i
139.         wm_r = COS(-2 * pi / m)
140.         wm_i = SIN(-2 * pi / m)
141.  
142.         FOR j = 0 TO n / 2 - 1 STEP m
143.             w_r = 1
144.             w_i = 0
145.  
146.             FOR k = 0 TO m / 2 - 1
147.                 p = j + k
148.                 q = p + (m \ 2)
149.  
150.                 u_r = w_r * xx_r(q) - w_i * xx_i(q)
151.                 u_i = w_r * xx_i(q) + w_i * xx_r(q)
152.                 v_r = xx_r(p)
153.                 v_i = xx_i(p)
154.  
155.                 xx_r(p) = v_r + u_r
156.                 xx_i(p) = v_i + u_i
157.                 xx_r(q) = v_r - u_r
158.                 xx_i(q) = v_i - u_i
159.  
160.                 u_r = w_r
161.                 u_i = w_i
162.                 w_r = u_r * wm_r - u_i * wm_i
163.                 w_i = u_r * wm_i + u_i * wm_r
164.             NEXT
165.         NEXT
166.     NEXT
167.  
168.     xx_r(n / 2) = xx_r(0)
169.     xx_i(n / 2) = xx_i(0)
170.  
171.     FOR i = 1 TO n / 2 - 1
172.         xx_r(n / 2 + i) = xx_r(n / 2 - i)
173.         xx_i(n / 2 + i) = xx_i(n / 2 - i)
174.     NEXT
175.  
176.     DIM xpr AS DOUBLE, xpi AS DOUBLE
177.     DIM xmr AS DOUBLE, xmi AS DOUBLE
178.  
179.     FOR i = 0 TO n / 2 - 1
180.         xpr = (xx_r(i) + xx_r(n / 2 + i)) / 2
181.         xpi = (xx_i(i) + xx_i(n / 2 + i)) / 2
182.  
183.         xmr = (xx_r(i) - xx_r(n / 2 + i)) / 2
184.         xmi = (xx_i(i) - xx_i(n / 2 + i)) / 2
185.  
186.         xx_r(i) = xpr + xpi * COS(2 * pi * i / n) - xmr * SIN(2 * pi * i / n)
187.         xx_i(i) = xmi - xpi * SIN(2 * pi * i / n) - xmr * COS(2 * pi * i / n)
188.     NEXT
189.  
190.     'symmetry, complex conj
191.     FOR i = 0 TO n / 2 - 1
192.         xx_r(n / 2 + i) = xx_r(n / 2 - 1 - i)
193.         xx_i(n / 2 + i) = -xx_i(n / 2 - 1 - i)
194.     NEXT
195. END SUB
196.  

[_vince]

Nice, looks interesting.  I'll have to update QB64 at some point to get the new memsound.  I wish I would get more involved in this but I haven't had much interest in qb64 recently.  Btw, I did not invent these algorithms, it's decades old math translated into qb64.

[Sanmayce]

Nice, looks interesting.  I'll have to update QB64 at some point to get the new memsound.  I wish I would get more involved in this but I haven't had much interest in qb64 recently.  Btw, I did not invent these algorithms, it's decades old math translated into qb64.

Hi Vince and Petr,
thanks for sharing this exciting etude.

Wanna write (with help and feedback from the community) a decent (at least not crappy) MP3/WAV frequency visualizer, will share my code in next days, at the moment I have a draft. The main goal is to have the jumping bars corresponding to the actual frequency peaks, to be not only synchronized but to drop no "frames".

By the way, @_vince , what was the source you used, can you share the link.
If your code fails to deliver, my intention is to write my own DFT etude, not FFT though - just the simplicity itself.

[Sanmayce]

Love seeing those jumping bars, I am eager to announce what is on my mind - e-pi-ALTES - the QB64 Audio_(Frequency_Spectrum)_Visualizer:

 

Three things inspired me to name my bar-jumper after the famous traitor:
- To honor Euler with 'e', to honor 'π' as well, as they are the basics of DFT (Discrete-Fourier-Transform);
- The etymology being "one who leaps upon";
- The superb dialogue in '300' movie, emphasizing the "beating heart", where the heart is the audio and spear/shield/sword are different frequencies:

Ephialtes : My father trained me to feel no fear to make spear and shield and sword as much a part of me as my own beating heart!
King Leonidas : [Ephialtes shows King Leonidas his thrust; it's good and the King is surprisingly impressed]  A fine thrust.

Indeed those jumping bars resemble the thrust of Epialtes' spear.

Note: Some write 'ph' and pronounce it as 'f', however the Greek historian Herodotus spelled it as Ἐπιάλτης, if you listen to the movie, the traitor himself pronounces P not F.

[_vince]

By the way, @_vince , what was the source you used, can you share the link.
If your code fails to deliver, my intention is to write my own DFT etude, not FFT though - just the simplicity itself.

https://qb64forum.alephc.xyz/index.php?topic=1938.msg111661#msg111661

[STxAxTIC]

I'll be keeping an eye on this thread. A while ago I almost had an oscilloscope working, but I think it suffered a windowing problem. Makes me wanna try it over again.... Good luck!

[MasterGy]

This is very good Petr! This code would have been good for me 15 years ago :)

[Petr]

Thanks, @MasterGy, but I've only slightly modified it. @_vince has the main merit.

[Sanmayce]

Glad to share the first DRAFT revision,
QB64 source plus the 32bit binary are attached.

 

The jumping of bars is too fast (even at 30 Frames-Per-Second), but this is RAW as it gets, for smoothing the animation some tweaks are needed, no idea at the moment what.

 

Critical feedback is welcome, also suggestions what can be bettered.

[Petr]

Hi, thanks for sharing. I just look in it (I have to drive with the presents for parrents in a while) and the analyzer, in order for you to go smoothly, the problem is here: (line around 359, I rummaged in it a bit):
         
                                              peakmagMatrix (Analyze) = peakmagMatrix (Analyze) * MultiplierForThePOWER

This is mathematically correct, but for the smooth running of the animation, I think I solved it so that if the old value is greater than the new one, then the old value will decrease, the correct low one will not be assigned immediately. I think that's how I dealt with it. I will still look at it, but only after Christmas.

I no longer remember whether I did the signal drop according to the new sample linearly or using percentages (when the difference between the old and new signal was for example 30 percent (power), X was then shifted 15 percent lower), but I don't know if I used this method in this matter.

[Sanmayce]

Glad to share the first operational revision, the QB64 source and the Windows executable (64bit, -O3 -mavx2) are in the attached file.
In order to play a WAV/MP3, just give the name as a parameter on the command line. Use quotes if necessary.

 

My i5-7200U struggled bigtime when changing the FFT samples from 1024 to 4096, at 8192 the tearing is nasty.
I have an idea to precalculate (in next revisions) all the 8192 samples long chunk and to JUST draw the magnitudes.

Also, I encountered a nasty noise-like problem in the decoded stream (bug?!), many sequences with unnatural high values making the visual output blinking badly, what is this, anyone?!

EDIT, 2021-Dec-29:
The bug was fixed, now the revision 2 is fully functional. I hate when don't have enough time to finish something, anyway, found it, the problem was in reading the decompressed stream with at offsets not divisible by 2, duh.

I tested e-pi-ALTES with Japanese narration (Sekirei, Ninja Resurrection) - the visualization is as it should - the vocal range is well presented, some nice bell shapes in the Cyan zone the 1KHz-5KHz, to me visualizer of the magnitudes of frequencies lacking in showing a RICH presence in vocal range is of no use, at least to me.

 

Visualizing the male vocal of the demon advisor Mori:

 

You may press Alt+Enter to toggle fullscreen/windowed modes, on my i5-7200U the transforming speed is at 25 frames-per-second for the highest FFT chunk - 8192 samples, it can vary between 128, 256, 512, 1024, 2048, 4096, 8192:

 

As always, the source code is in the attached package, also 32bit and 64bit executables for Windows.

Enjoy!

[Sanmayce]

Glad to share the latest and finest e-pi-ALTES, it features the "snowfall" effect:

 

The rewind (backwards or forwards) can be made with Left/Right arrow as well.
The volume Up/Down can be made with Up/Down arrow as well.

Also, did some screen-resolution-aware actions, on old laptops (with 1366x768) e-pi-ALTES enforces the 8x8 Toshiba font, on FHD (1920x1080) laptops it enforces another superb Japanese font but 6x12, my calculations are in the source code as comments as to how to preserve undistorted scaling when you go FULL-SCREEN, hope it looks good.

I strongly recommend using e-pi-ALTES when listening to audiobooks or Japanese audio narration, the later are AWESOMELY done, visualization comes ... falling, as in the superb song ending the classic film "Roadhouse", it goes "The night comes falling".

I love this very revision, for the first time I see what I hear, hee-hee.

Having compiled it on Linux Fedora 35, here comes the promo video (this is how old laptops will look like, in 1680 or 1366 modes):

Have a nice new year!