PureBasic

Forums PureBasic - Français
http://forums.purebasic.com/french/

Gestion des fréquences

http://forums.purebasic.com/french/viewtopic.php?t=4373

Page 1 sur 1

Gestion des fréquences

Publié : mer. 01/févr./2006 11:40
par Pascal Vernie
Bonjour à tous
Je suis en train de faire un lecteur MP3 cela va pas trop mal pour l'instant.
Mais je voudrais lui ajouter pour l'esthetique des barres qui modulent en fonction de la fréquence (Par Ex 100Hz,250Hz,500Hz etc).
J'ai cherché sur le forum et sauf mauvaise recherche de ma part je n'ai rien trouvé, quelqu'un pourrait-il me donner quelques pistes sur le sujet.
En vous remerciant d'avance.
Pascal

Publié : mer. 01/févr./2006 12:17
par Progi1984
Essaie de voir dans la droopy lib, je crois qu'il ya quelque chose qui pourrait t'intéresser !

Publié : mer. 01/févr./2006 12:41
par KarLKoX
J'ai fait une routine FFT via la lib F64 (donc 64 bits), je ne l'ai pas essayé à fond mais pour un filtre passe bas, ça passe donc pour un simple spectre ça devrait le faire. Il n'y a aucune doc ni quoi que ce soit, il faudra que tu lise un peu le code pour l'utiliser (rien de bien compliquer en soit) mais en gros, il suffira juste que tu envoye un bout de donnée sonore à filtrer, ça te renverra le FFT de celui ci, il faudra que tu fasse un log (ou eq pour un truc tip top) de la somme des données / 2. (nyquist tout ça)
J'Up ça dès que je rentre du boulot (çe soir).

Publié : mer. 01/févr./2006 20:18
par KarLKoX

Code : Tout sélectionner

;- FFT : THE THEORY
;    This computes an in-place complex-To-complex FFT
;    x And y are the real And imaginary arrays of 2^m points.
;    dir =  1 gives forward transform
;    dir = -1 gives reverse transform
; 
;      Formula: forward
;                   N-1
;                   ---
;               1   \          - j k 2 pi n / N
;       X(n) = ---   >   x(k) e                    = forward transform
;               N   /                                n=0..N-1
;                   ---
;                   k=0
; 
;       Formula: reverse
;                   N-1
;                   ---
;                   \          j k 2 pi n / N
;       X(n) =       >   x(k) e                    = forward transform
;                   /                                n=0..N-1
;                   ---
;                   k=0

;- WRITE FFT VALUES TO A FILE
Procedure WriteSeries_F64(fname.s,*x.double,*y1.double,*y2.double,n.l)
Protected handle.l
Protected values.s
Protected tmpValue.double, floatValue.f

	handle = CreateFile(#PB_Any, fname)
	If handle = 0
		MessageRequester("Error", "An error occured when creating the file " + fname)
		ProcedureReturn #False
	EndIf
	
	For i = 0 To n - 1
		CopyMemory(*x+i*SizeOf(double), tmpValue, SizeOf(double) )
		floatValue = F64_toFloat(tmpValue)	
		values = Str(i) + " " + StrF(floatValue) + " " 
		CopyMemory(*y1+i*SizeOf(double), tmpValue, SizeOf(double) )		
		floatValue = F64_toFloat(tmpValue)		
		values + StrF(floatValue) + " " 
		CopyMemory(*y2+i*SizeOf(double), tmpValue, SizeOf(double) )		
		floatValue = F64_toFloat(tmpValue)		
		values + StrF(floatValue) + Chr(13)
						
		WriteString(values)
	Next i
	CloseFile(handle)
	
	ProcedureReturn #True		
EndProcedure

;- FFT FUNCTION
Procedure FFT_F64(dir.l, m.l, *x.double, *y.double)
Protected nn, i, i1, j, k, i2, l1, l2
Protected c1.double, c2.double, tx.double, ty.double, t1.double, t2.double, u1.double, u2.double, z.double
Protected tmpValue1.double, tmpValue2.double
Protected tmpValue3.double, tmpValue4.double
Protected tmpValue5.double

	; Calculate the number of points
	nn = 1
	For i=0 To m - 1
      nn * 2
	Next i

	; Do the bit reversal
	i2 = nn >> 1
	j = 0
	For i=0 To nn - 1 - 1
		If i < j
 	  	CopyMemory(*x+i*SizeOf(double), tx, SizeOf(double) )									; tx 	= x[i]
   		CopyMemory(*y+i*SizeOf(double), ty, SizeOf(double) )									; ty 	= y[i]
 		CopyMemory(*x+j*SizeOf(double), *x+i*SizeOf(double), SizeOf(double) )		; x[i] = x[j]
 		CopyMemory(*y+j*SizeOf(double), *y+i*SizeOf(double), SizeOf(double) )		; y[i] = y[j]
   		CopyMemory(tx, *x+j*SizeOf(double), SizeOf(double) )									; x[j] = tx
 	  	CopyMemory(ty, *y+j*SizeOf(double), SizeOf(double) )									; y[j] = ty
		EndIf
      	k = i2
		While k <= j
			j - k
			k >> 1
		Wend
      	j + k
	Next i

	; Compute the FFT
	F64_Float(c1, -1.0)
	F64_Float(c2, 0.0)
	l2 	= 1
	For l=0 To m - 1
		l1 	= l2
		l2 	<< 1
   		F64_Float(u1, 1.0)
   		F64_Float(u2, 0.0)
		F64_Float(tmpValue1, 0.0)
  		F64_Float(tmpValue2, 0.0)
		F64_Float(tmpValue3, 0.0)
  		F64_Float(tmpValue4, 0.0)
		F64_Float(tmpValue5, 0.0)
		For j = 0 To l1 - 1
         	For i = j To nn - 1
				    i1 = i + l1
				    CopyMemory(*x+i1*SizeOf(double), tmpValue1, SizeOf(double) )    	; tmpValue1 = x[i1]
				    F64_Mul(u1, tmpValue1, tmpValue3)											 	; tmpValue3 = u1 * x[i1]
				    CopyMemory(*y+i1*SizeOf(double), tmpValue2, SizeOf(double) )	 	; tmpValue2 = y[i1]
				    F64_Mul(u2, tmpValue2, tmpValue4)											 	; tmpValue4 = u2 * y[i1]				
				    F64_Sub(tmpValue3, tmpValue4, t1)											 	; t1 = u1 - u2
				
				    F64_Mul(u2, tmpValue1, tmpValue3)											 	; tmpValue3 = u2 * x[i1]
				    F64_Mul(u1, tmpValue2, tmpValue4)												; tmpValue4 = u1 *[y1]
				    F64_Add(tmpValue3, tmpValue4, t2)												; t2 = u1 + u2
				
				    CopyMemory(*x+i*SizeOf(double), tmpValue3, SizeOf(double) )		  	; tmpValue1 = x[i]
				    F64_Sub(tmpValue3, t1, tmpValue1)												; tmpValue1 = tmpValue3 - t1
				    CopyMemory(tmpValue1, *x+i1*SizeOf(double), SizeOf(double) )	  	; x[i1] = tmpValue1
    				
				    CopyMemory(*y+i*SizeOf(double), tmpValue4, SizeOf(double) )		  	; tmpValue2 = y[i]
				    F64_Sub(tmpValue4, t2, tmpValue1)												; tmpValue1 = tmpValue4 - t2
				    CopyMemory(tmpValue1, *y+i1*SizeOf(double), SizeOf(double) )	  	; y[i1] = tmpValue4
								
				    F64_Add(tmpValue3, t1)													                ; tmpValue3 = tmpValue3 + t1
				    CopyMemory(tmpValue3, *x+i*SizeOf(double), SizeOf(double) )	    	; x[i] = tmpValue3
				
				    F64_Add(tmpValue4, t2)  	  											              	; tmpValue4 = tmpValue4 + t2
				    CopyMemory(tmpValue4, *y+i*SizeOf(double), SizeOf(double) )	    	; y[i] = tmpValue4				        
	          	          
          	i + l2 - 1	; step l2
	          
         	Next i
       	         																							  	 	    ; (x1 and x2 are here to simplify what i am doing)
			  F64_Mul(u1, c1, tmpValue1)																	; x1 = u1 * c1
			  F64_Mul(u2, c2, tmpValue2)															  		; x2 = u2 * c2				
			  F64_Sub(tmpValue1, tmpValue2, z)														; z = x1 - x2
			  F64_Mul(u1, c2, tmpValue1)																  	; x1 = u1 * c2
			  F64_Mul(u2, c1, tmpValue2)																	; x2 = u2 * c1
			  F64_Add(tmpValue1, tmpValue2, u2)													; u2 = x1 + x2
			  F64_Copy(u1, z)									  												; u1 = z
  			
		Next j
		F64_Float(tmpValue1, 1.0)																			; tmpValue1 = 1.0
		F64_Float(tmpValue2, 2.0)	  	 																	; tmpValue2 = 2.0		
		F64_Sub(tmpValue1, c1, tmpValue3)  														; tmpValue3 = 1.0 - c1
		F64_Div(tmpValue3, tmpValue2)																	; tmpValue3 =  tmpValue3 / 2.0
		F64_Sqrt(tmpValue3, c2)																			; c2 = sqrt(tmpValue3)
      	If (dir = 1)
        	 F64_Negate(c2)                                                   									; c2 = -c2
      	EndIf
		F64_Add(tmpValue1, c1, tmpValue3)															; tmpValue3 = 1.0 + c1
		F64_Div(tmpValue3, tmpValue2)																	; tmpValue3 = tmpValue3 / 2.0
		F64_Sqrt(tmpValue3, c1)					  														; c1 = sqrt(tmpValue3)	      
		
	Next l

	; Scaling For forward transform
	If dir = 1
		For i = 0 To nn - 1
			CopyMemory(*x+i*SizeOf(double), tmpValue1, SizeOf(double) )				    ; tmpValue1 = x[i]	
			F64_Int(tmpValue2, nn)																			; tmpValue2 = nn
			F64_Div(tmpValue1, tmpValue2)																; tmpValue1 / nn
			CopyMemory(tmpValue1, *x+i*SizeOf(double), SizeOf(double) )				    ; x[i] = tmpValue1
			
			CopyMemory(*y+i*SizeOf(double), tmpValue1, SizeOf(double) )				    ; tmpValue1 = y[i]
			;*tmpValue5 = *y+i*SizeOf(double)                                     						; tmpValue5 = y[i]
			F64_Div(tmpValue1, tmpValue2)																; tmpValue2 / nn
			CopyMemory(tmpValue1, *y+i*SizeOf(double), SizeOf(double) )				    ; y[i] = tmpValue1			
		Next i
	EndIf

   ProcedureReturn #True

EndProcedure

;- FFT EXAMPLE

; Test-bed For Fourier digital filter design
; Illustrative code only....not a "product".


#N					=	1024					; Number of samples
#NP			  		= 10					  	; Power of 2, 2^NP = N
#FREQ		  		=	500.0				; Sampling frequency
#F1					=	50.0					; Cutoff frequency
#W					=	20					; Desired filter width
#RECT		  		=	0
#HANN		  		=	1
#WINDOW   		=	1

#TWOPI			= 6.28318530

#HEAP_GENERATE_EXCEPTIONS 		= 4
#HEAP_ZERO_MEMORY					  	= 8

Procedure AllocOK(*ptr)
	If *ptr = 0
		MessageRequester("Error", "Memory allocation failed")
		End
	EndIf
EndProcedure

Procedure DoFFT_F64()
Protected i
Protected *re.double, *im.double, *freq.double, *time.double, *window.double, *mag.double, *pha.double
Protected tmpValue1.double, tmpValue2.double
Protected tmpValue3.double, tmpValue4.double
Protected tmpValue5.double, tmpValue6.double

  *re = AllocateMemory(#N * SizeOf(double) )
  AllocOK(*re)
  *im = AllocateMemory(#N * SizeOf(double) )
  AllocOK(*im)
  *freq = AllocateMemory(#N * SizeOf(double) )
  AllocOK(*freq)
  *time = AllocateMemory(#N * SizeOf(double) )
  AllocOK(*time)
  *window = AllocateMemory(#N * SizeOf(double) )
  AllocOK(*window)
  *mag = AllocateMemory(#N * SizeOf(double) )
  AllocOK(*mag)
  *pha = AllocateMemory(#N * SizeOf(double) )
  AllocOK(*pha)
  
	; Fill in the time And frequency indices
	; These are only used To give more informative file output
	; although they are useful To ensure you understand the Data arrangement
	For i = 0 To #N - 1	
		F64_Float(tmpValue1, i)
		F64_Int(tmpValue2, #FREQ)
		F64_Div(tmpValue1, tmpValue2)
		CopyMemory(tmpValue1, *time+i*SizeOf(double), SizeOf(double) )		
		
		F64_Float(tmpValue1, #N)
		F64_Val(tmpValue2, "2.0")
		F64_Div(tmpValue1, tmpValue2)
		F64_Int(tmpValue3, i)
		If F64_Cmp(tmpValue3, tmpValue1) = -1 Or F64_Cmp(tmpValue3, tmpValue1) = 0	; Positive frequencies
			F64_Int(tmpValue1, #FREQ)			
			F64_Int(tmpValue2, #N)
			F64_Div(tmpValue1, tmpValue2)
			F64_Mul(tmpValue1, tmpValue3)
			CopyMemory(tmpValue1, *freq+i*SizeOf(double), SizeOf(double) )	
		Else															                                          ; Negative frequencies
			F64_Int(tmpValue1, i)			
			F64_Int(tmpValue2, #FREQ)
			F64_Mul(tmpValue1, tmpValue2)
			F64_Float(tmpValue3, #N)
			F64_Div(tmpValue1, tmpValue3)
			F64_Sub(tmpValue1, tmpValue2)
			CopyMemory(tmpValue1, *freq+i*SizeOf(double), SizeOf(double) )				
		EndIf			
	Next i

	; Create the impulse response of the filter in the frequency domain 
	; This is an example of a low pass filter, cutoff at F1 Hz
	; Note that since I am using a full complex FFT it is necessary
	; To correctly fill in the negative frequencies as well.
	; Zero phase is assumed at this stage but this will become linear
	; phase when the filter is shifted in time to make it causal.
	For i = 0 To #N - 1
		F64_Float(tmpValue1, 0.0)
		CopyMemory(tmpValue1, *re+i*SizeOf(double), SizeOf(double) )		
 		;If i <= (#F1 * #N / #FREQ) Or i >= (#N - #F1 * #N / #FREQ)
		F64_Float(tmpValue1, #F1)
		F64_Float(tmpValue2, #N)
		F64_Mul(tmpValue1, tmpValue2)
		F64_Float(tmpValue2, #FREQ)
		F64_Div(tmpvalue1, tmpValue2)
		
		F64_Float(tmpValue3, #F1)
		F64_Float(tmpValue4, #N)
		F64_Mul(tmpValue3, tmpValue4)
		F64_Float(tmpValue4, #FREQ)
		F64_Div(tmpValue3, tmpValue4)
		F64_Float(tmpValue4, #N)
		F64_Sub(tmpValue4, tmpValue3, tmpValue3)
		F64_Int(tmpValue2, i)
		If F64_Cmp(tmpValue2, tmpValue1) = - 1 Or F64_Cmp(tmpValue2, tmpValue3) = 1
			F64_Float(tmpValue2, 1.0)
			CopyMemory(tmpValue2, *re+i*SizeOf(double), SizeOf(double) )
 		EndIf			
		F64_Float(tmpValue2, 0.0)
		CopyMemory(tmpValue2, *im+i*SizeOf(double), SizeOf(double) )
	Next i
	WriteSeries_F64("impulse1.f",*freq,*re,*im,#N)
			
 	; Inverse FFT the signal into the time domain
 	; If the filter is created properly then the imaginary
 	; components should be zero....baring numerical error.
	FFT_F64(-1,#NP,*re,*im)
	WriteSeries_F64("impulse.t",*time,*re,*im,#N)
 		
  ; Create the window so to make the series finite duration
	For i = 0 To #N - 1
    F64_Float(tmpValue1, 0.0)
		CopyMemory(tmpValue1, *window+i*SizeOf(double), SizeOf(double) )		
		F64_Int(tmpValue1, i)
		F64_Int(tmpValue2, #W)
		F64_Int(tmpValue3, 2)
		F64_Int(tmpValue4, #N)
		F64_Div(tmpValue2, tmpValue3, tmpValue5) ; tmpValue5 = #W / 2
		;If i <= #W / 2
		If F64_Cmp(tmpValue1, tmpValue5) = -1 Or F64_Cmp(tmpValue1, tmpValue5) = 0
			If (#WINDOW = #RECT)
				F64_Float(tmpValue6, 1.0)
				CopyMemory(tmpValue6, *window+i*SizeOf(double), SizeOf(double) )	
			Else
				F64_Float(tmpValue6, 0.5)
				F64_PI(tmpValue4)
				F64_Mul(tmpValue4, tmpValue3) ; tmpValue4 = #TWOPI
				F64_Mul(tmpValue1, tmpValue4) ; tmpValue1 = i * #TWOPI
				F64_Div(tmpValue1, tmpValue2) ; tmpValue1 / W
				F64_Cos(tmpValue1)            ; tmpValue1 = Cos(tmpValue1)
				F64_Mul(tmpValue1, tmpValue6) ; tmpValue1 * 0.5
				F64_Add(tmpValue1, tmpValue6) ; tmpValue1 = tmpValue1 + 0.5
				CopyMemory(tmpValue1, *window+i*SizeOf(double), SizeOf(double) )	
			EndIf				
		F64_Sub(tmpValue4, tmpValue5, tmpValue6)
		;ElseIf (i >= #N - #W/2) 		
		ElseIf F64_Cmp(tmpValue1, tmpValue6) = 1 Or F64_Cmp(tmpValue1, tmpValue6) = 0
      If (#WINDOW = #RECT)
        F64_Float(tmpValue6, 1.0)
        CopyMemory(tmpValue6, *window+i*SizeOf(double), SizeOf(double) )	
      Else
        F64_Float(tmpValue6, 0.5)
				F64_PI(tmpValue5)
				F64_Mul(tmpValue5, tmpValue3) ; tmpValue5 = #TWOPI
				F64_Div(tmpValue5, tmpValue2) ; tmpValue5 / #W
  			F64_Sub(tmpValue4, tmpValue1) ; tmpValue4 = #N - i
				F64_Mul(tmpValue4, tmpValue5) ; tmpValue4 = (#N - i) * #TWOPI / #W
				F64_Cos(tmpValue4)            ; tmpValue4 = Cos(tmpValue4)
  			F64_Mul(tmpValue4, tmpValue6) ; tmpValue4 * 0.5
				F64_Add(tmpValue4, tmpValue6) ; tmpValue4 = tmpValue4 + 0.5
  			CopyMemory(tmpValue4, *window+i*SizeOf(double), SizeOf(double) )	
			EndIf				
		EndIf
	Next i
	WriteSeries_F64("window.t",*time,*window,*im,#N)
	
  ; Apply the window
	For i = 0 To #N - 1
	  CopyMemory(*window+i*SizeOf(double), tmpValue1, SizeOf(double) )
  	CopyMemory(*re+i*SizeOf(double), tmpValue2, SizeOf(double) )
	  CopyMemory(*im+i*SizeOf(double), tmpValue3, SizeOf(double) )
		F64_Mul(tmpValue1, tmpValue2, tmpValue4)
		F64_Float(tmpValue5, 0.0)
		CopyMemory(tmpValue4, *re+i*SizeOf(double), SizeOf(double) )
		CopyMemory(tmpValue5, *im+i*SizeOf(double), SizeOf(double) )
	Next i
	WriteSeries_F64("windowed.t",*time,*re,*im,#N)
		
	; Shift the filter to make it causal
	For i = #W To 0 Step - 1
		CopyMemory(*re+i*SizeOf(double), tmpValue1, SizeOf(double) )
		re_index.l = i - #W / 2		
		CopyMemory(*re+re_index*SizeOf(double), tmpValue2, SizeOf(double) )
		CopyMemory(tmpValue2, *re+i*SizeOf(double), SizeOf(double) )
	Next i		
	For i = 0 To #W / 2 - 1
		CopyMemory(*re+i*SizeOf(double), tmpValue1, SizeOf(double) )
		re_index.l = #N - #W/2 + i
		CopyMemory(*re+re_index*SizeOf(double), tmpValue2, SizeOf(double) )
		CopyMemory(tmpValue2, *re+i*SizeOf(double), SizeOf(double) )
		F64_Float(tmpValue3, 0.0)
		CopyMemory(tmpValue3, *re+re_index*SizeOf(double), SizeOf(double) )
	Next i
	WriteSeries_F64("causal.t",*time,*re,*im,#N)
	
	; Transform back into the frequency domain mostly as a check
	FFT_F64(1,#NP,*re,*im)
	For i = 0 To #N - 1
		CopyMemory(*re+i*SizeOf(double), tmpValue1, SizeOf(double) )
		CopyMemory(*im+i*SizeOf(double), tmpValue2, SizeOf(double) )
		F64_Mul(tmpValue1, tmpValue1, tmpValue3)  ; tmpValue3 = re[i] * re[i]
		F64_Mul(tmpValue2, tmpValue2, tmpValue4)  ; tmpValue4 = im[i] * im[i]
		F64_Add(tmpValue3, tmpValue4, tmpValue5)  ; tmpValue3 = re[i]*re[i] + im[i]*im[i]
		F64_Sqrt(tmpValue5)
		CopyMemory(tmpValue5, *mag+i*SizeOf(double), SizeOf(double) )
		
		F64_AtanQuad(tmpValue2, tmpValue1, tmpValue3)
		CopyMemory(tmpValue3, *pha+i*SizeOf(double), SizeOf(double) )
	Next i
	WriteSeries_F64("impulse2.f",*freq,*re,*im,#N)
	WriteSeries_F64("impulse3.f",*freq,*mag,*pha,#N)
	
	FreeMemory(*re)
	FreeMemory(*im)
	FreeMemory(*freq)
	FreeMemory(*time)
	FreeMemory(*window)
	FreeMemory(*mag)
	FreeMemory(*pha)
		
EndProcedure

DoFFT_F64()
Bon, l'exemple est bidon (il est même pas bon) mais la routine FFT elle, ne l'est pas et fonctionne tres bien.
Tu pourras trouver facilement du code sur le net pour utiliser les FFT pour un petit eq.
Bonne chance ;)

Publié : mer. 01/févr./2006 20:49
par Flype
Je crois que tu vas bientôt être obligé ( ou tenté ) de recoder tout çà sans la lib F64, avec PB4.0 :D
C'est intéressant ton truc, ce serait plutot cool d'inclure çà dans le recorder de SoundEditor (http://www.freesoundeditor.com/incagen.html).

Publié : mer. 01/févr./2006 21:08
par KarLKoX
Oui, j'ai d'autres routines FFT en cours (via les pointeurs, fonctions pb native et un embryon de simd).
En fait, tout ce que je fait en ce moment est pensé pour être implémenté facilement dans soundeditor donc dès que j'aurais réussi à faire ce que je souhaite, j'ajouterais ça petit à petit :)

Publié : mer. 01/févr./2006 21:23
par Flype
Cool, je t'encourage. :twisted:

Publié : mer. 01/févr./2006 22:12
par Droopy
Comme l'indique Progi, pour créer les barres ( vu-mêtre ), dans la Droopy Lib les fonctions Bargraph devraient t'aider.

Publié : mer. 01/févr./2006 23:35
par Pascal Vernie
Bonsoir
Je vous remercie pour vos réponses.
J'ai vu effectivement dans la lib de Droopy pour l'affichage cela me semble parfait pour ce que je veux faire.
Pour KarLKoX j'ai du pain sur la planche. Je suis pas encore un as du PureBasic, mais avec de la persévérance on devrait y arriver...........
Merci encore de votre aide.
Pascal
Fuseau horaire sur UTC+01:00
Page 1 sur 1

Développé par phpBB® Forum Software © phpBB Limited

Traduction française officielle © Qiaeru