Gestion des fréquences

Pascal Vernie · Message par **Pascal Vernie** » mer. 01/févr./2006 11:40

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

Progi1984 · Message par **Progi1984** » mer. 01/févr./2006 12:17

Essaie de voir dans la droopy lib, je crois qu'il ya quelque chose qui pourrait t'intéresser !

KarLKoX · Message par **KarLKoX** » mer. 01/févr./2006 12:41

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).

KarLKoX · Message par **KarLKoX** » mer. 01/févr./2006 20:18

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

Flype · Message par **Flype** » mer. 01/févr./2006 20:49

Je crois que tu vas bientôt être obligé ( ou tenté ) de recoder tout çà sans la lib F64, avec PB4.0

C'est intéressant ton truc, ce serait plutot cool d'inclure çà dans le recorder de SoundEditor (http://www.freesoundeditor.com/incagen.html).

KarLKoX · Message par **KarLKoX** » mer. 01/févr./2006 21:08

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

Flype · Message par **Flype** » mer. 01/févr./2006 21:23

Cool, je t'encourage.

Droopy · Message par **Droopy** » mer. 01/févr./2006 22:12

Comme l'indique Progi, pour créer les barres ( vu-mêtre ), dans la Droopy Lib les fonctions Bargraph devraient t'aider.

Pascal Vernie · Message par **Pascal Vernie** » mer. 01/févr./2006 23:35

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