Visual Basic 6.0 - Preview:
'====================================================================== ' Descrizione.....: Collezione di routines e costanti di utilita' ' per il Digital Signal Processing. '====================================================================== ' ' Le routines di questo modulo sono un sotto insieme ' di quelle contenute nel progetto FiltCalc (sito dei ' DownLoads) e sono state modificate per adattarle a ' questa applicazione. ' Non tutte le costanti e le routines di questo ' modulo vengono, necessariamente, usate. ' Option Explicit ' '--- Windows ---------------------------------------------------------- Public Type Window_Type ' Nome As String ' Nome della "Window": PMin As Double ' Valore Min. del Parametro associato. PMax As Double ' Valore Max. del Parametro associato. PCor As Double ' Valore corrente del Parametro associato. End Type ' Private Const A0# = 0.99938 ' Coefficienti per la Private Const A1# = 0.041186 ' Weber Window. Private Const A2# = -1.637363 ' Private Const A3# = 0.828217 ' Private Const B0# = 1.496611 ' Private Const B1# = -1.701521 ' Private Const B2# = 0.372793 ' Private Const B3# = 0.0650621 ' ' '--- Spettro mobile (sliding spectrum) ------------------------------------------------- ' Implementazione dei circuiti tratti da: ' "Theory and Application of Digital Signal Processing" ' di L. Rabiner e B. Gold. - pg. 382, 383. ' Private NsSpm& ' N° di campioni per il calcolo dello spettro mobile. Private NormSpm# ' Fattore di normalizzazione sul N° di campioni. ' ' Coefficienti per il calcolo con convoluzione diretta: Private z1() As Complex ' vettore dei coefficienti di calcolo dello spettro mobile. Private SmRE#() ' registro a scorrimento dei campioni del segnale. ' ' Coefficienti per il calcolo ricorsivo: Private Sn_1 As Complex ' Sn * Z^-1. Private z1_1 As Complex ' z1(0)^-1 Private z1_N As Complex ' z1(0)^-N Private Ic&, Kc& ' indici in SmRE() per buffer circolare. Private NsSpm1& ' NsSpm + 1. Private Xz1_N As Complex ' variabili di appoggio. Private Yz1_1 As Complex ' " " " Private SnI As Complex ' " " " Private Sn As Complex ' " " " ' ' Per multi frequenza: Private NFreqSPM_1& ' N° di frequenze - 1 a cui calcolare gli spettri. ' Coefficienti per il calcolo ricorsivo multi frequenza: Private Sn_MF_1() As Complex ' vettore degli (NFreqSPM_1 + 1) Sn * Z^-1 Private z1_MF_1() As Complex ' z1()^-1 Private z1_MF_N() As Complex ' z1()^-N ' Private SqRE#() ' Registro a scorrimento dei campioni del segnale al quadrato. Private VMed# ' Valore medio corrente. ' '--- Filtri FIR ------------------------------------------------------------------------ ' Algoritmi tradotti ed adattati dal FORTRAN di: ' "Digital Filters and their Applications" ' di V. Cappellini, A. G. Constantinides, P. Emiliani. ' Window method: pg. 350. ' Private Hc!() ' Vettore dei coefficienti del filtro. ' Private WF#() ' Tabella dei coefficienti per Windowing. ' Private SRE!() ' Registro a scorrimento dei Dati da Filtrare. ' Private Const FiltroErrFIR$ = "Le Routines di calcolo hanno trovato" & vbNewLine _ & "una condizione imprevista." & vbNewLine _ & "Provare a rivedere i parametri del filtro." & vbNewLine & vbNewLine ' '--- Filtri IIR ------------------------------------------------------------------------ ' Gli algoritmi per la sintesi di filtri del tipo Butterworth ' e Chebyshev sono stati tradotti ed adattati dal FORTRAN di: ' "Digital Filters and their Applications" ' di V. Cappellini, A. G. Constantinides, P. Emiliani. ' pg. 367, 368, 369. ' Metodi di calcolo dei filtri elementari. ' - Algoritmi di trasformazione da "Digital Signal Processing" ' di W. D. Stanley - pg. 172, 173, 174. ' - La funzione di trasferimento del risuonatore reale parallelo ' e' tratta da "Teoria delle Reti Elettriche", Appunti dai ' corsi del Politecnico di Torino - pg. (1.3)1 e seg. ' - L' idea della sostituzione degli zeri per il "Notch Filter" ' proviene da: http://www-users.cs.york.ac.uk/~fisher/mkfilter/res.html ' Private NK&, NCel& ' Ordine e numero di sezioni del filtro. Private Ac!() ' Coefficienti del filtro. Private Bc!() ' " " " ' Private Const NCMax& = 20 ' N. Massimo di sezioni del filtro. ¦ ' Private w!() ' Registri delle sezioni del filtro. ' Private CEB#(1 To 2 * NCMax + 1) ' Vettori in uso durante la sintesi. Private AN#(0 To 4, 0 To 2) ' " " " " Private FINA#(0 To 2) ' " " " " Private FINB#(0 To 2) ' " " " " ' Private Const FiltroErrIIR$ = "Le Routines di calcolo hanno trovato" & vbNewLine _ & "una condizione imprevista." & vbNewLine _ & "Provare a rivedere i parametri del filtro." & vbNewLine & vbNewLine ' '--- Per routines SFFTBI e SFTTBF ------------------------------------------------------ Private MM& ' integer such that N = 2**MM Private S1#() ' array of sin() table (length >= n/8-1) Private C1#() ' array of cos() table (length >= n/8-1) Private S3#() ' array of sin() table (length >= n/8-1) Private C3#() ' array of cos() table (length >= n/8-1) Private ITAB&() ' integer bit reversal table (length >= sqrt(2n)) Private D1#() ' Vettore dei dati di ingresso a base 1, come richiesto dal FORTRAN. ' '--- Per auto e mutua correlazione con SFFTBI_Corr, SFTTBF_Corr e SFFTBB_Corr ---------- Private N1_C& ' N° di valori in D1(). Private N2_C& ' N° di valori in D2(). Private MM_C& ' MM_C = Ceil(Log(CDbl(N1_C + N2_C - 1)) / Log2) Private N_C& ' N_C = 2**MM_C Private S1_C#() ' array of sin() table (length >= N_C/8-1) Private C1_C#() ' array of cos() table (length >= N_C/8-1) Private S3_C#() ' array of sin() table (length >= N_C/8-1) Private C3_C#() ' array of cos() table (length >= N_C/8-1) Private ITAB_C&() ' integer bit reversal table (length >= sqrt(2*N_C)) Private D1_C#() ' Vettore dei dati di ingresso a base 1, come richiesto dal FORTRAN. Private NFre_C& ' N° di frequenze calcolate. Private WnRe#() ' Tavole dei seni/coseni per il calcolo della trasformata del Private WnIm#() ' segnale ritardato di N1_C campioni (solo per AutoCorr_FT). Private Re1_C#() ' Parte reale della 1° trasformata. Private Im1_C#() ' " immaginaria " " " Private Re2_C#() ' Parte reale della 2° trasformata. Private Im2_C#() ' " immaginaria " " " Private s#(), f#() ' Vettori d' appoggio. ' '--- API di gestione memoria: ---------------------------------------------------------- Private Declare Sub CopyMemory Lib "kernel32" Alias "RtlMoveMemory" _ (hpvDest As Any, hpvSource As Any, ByVal numBytes As Long) ' Private Declare Sub MoveMemory Lib "kernel32" Alias "RtlMoveMemory" _ (hpvDest As Any, hpvSource As Any, ByVal numBytes As Long) ' Private Declare Sub ZeroMemory Lib "kernel32" Alias "RtlZeroMemory" _ (hpvDest As Any, ByVal numBytes As Long) Public Function EliminaCC(d() As Double, Optional ByRef VMed As Double) As Double() ' ' Ritorna i dati del vettore D(0 To N) dopo l' eliminazione ' della componente continua. I dati risultanti avranno, quindi, ' valor medio = 0. Ritorna anche VMed, valor medio dei dati in D(). ' Dim I&, N&, dd#() ' dd() = d() N = UBound(dd) VMed = 0# ' For I = 0 To N VMed = VMed + dd(I) Next I VMed = VMed / (N + 1) ' For I = 0 To N dd(I) = dd(I) - VMed Next I ' EliminaCC = dd() ' ' ' End Function Public Function SPM(ByVal Xn!) As Single ' ' Funzione per il calcolo con convoluzione diretta dello ' spettro mobile (sliding spectrum) su NsSpm campioni. ' Devono essere stati definiti, precedentemente, il vettore ' dei coefficienti Z1() ed un vettore SmRE(0 To NsSpm - 1) da ' usare come registro a scorrimento dei campioni del segnale. ' Dim K& ' For K = NsSpm - 1 To 1 Step -1 SmRE(K) = SmRE(K - 1) Next K SmRE(0) = CDbl(Xn) ' Sn = CCmp(SmRE(0), 0#) For K = 1 To NsSpm - 1 Sn = CSom(Sn, CMol(CCmp(SmRE(K), 0#), z1(K))) Next K ' SPM = CSng(CAbs(Sn)) / CDbl(NsSpm) ' ' ' End Function Public Sub SPM_Init(ByVal NsSpm_I As Long, ByVal Freq As Double, ByVal Fs As Double) ' ' Calcolo dei coefficienti necessari alla ' implementazione con convoluzione diretta ' dello spettro mobile su NsSpm campioni: ' Dim I&, OmegaT# ' NsSpm = NsSpm_I ' N° di campioni per il calcolo dello spettro mobile. NormSpm = CDbl(NsSpm) / 2# ' Fattore di normalizzazione sul N° di campioni. ' ReDim z1(0 To NsSpm) ' Vettore dei coefficienti di calcolo con ' convoluzione diretta dello spettro mobile; ' z1(0) = e^j*omega*T. ReDim SmRE#(0 To NsSpm) ' Registro a scorrimento dei campioni del segnale. ' OmegaT = PI2 * Freq / Fs ' Per frequenza in [Hz]. 'OmegaT = PI2 * Freq ' Per frequenza in [f/fs]. z1(0) = CExp(CCmp(0#, OmegaT)) ' ' Per calcolo con convoluzione diretta: For I = 1 To NsSpm - 1 z1(I) = CPtN(z1(0), -I) Next I ' ' ' End Sub Public Sub SPM_R_Init(ByVal NsSpm_I As Long, ByVal Freq As Double, ByVal Fs_I As Double) ' ' Calcolo dei coefficienti necessari alla ' implementazione ricorsiva dello spettro ' mobile su NsSpm campioni: ' Dim I&, OmegaT# ' NsSpm = NsSpm_I ' N° di campioni per il calcolo dello spettro mobile. NsSpm1 = NsSpm + 1 NormSpm = CDbl(NsSpm) / 2# ' Fattore di normalizzazione sul N° di campioni. ' ReDim z1(0 To 0) ' z1(0) = e^j*omega*T. ReDim SmRE#(0 To NsSpm) ' Registro a scorrimento dei campioni del segnale. ' OmegaT = PI2 * Freq / Fs_I ' Per frequenza in [Hz]. 'OmegaT = PI2 * Freq ' Per frequenza in [f/fs]. z1(0) = CExp(CCmp(0#, OmegaT)) ' ' Per calcolo ricorsivo: Sn_1 = CCmp(0#, 0#) z1_1 = CPtN(z1(0), -1) z1_N = CPtN(z1(0), -NsSpm) ' ' ' End Sub Public Sub SPM_MF_R_Init(ByVal NsSpm_I As Long, ByVal NFRE As Long, Freq() As Double, _ ByVal Fs As Double) ' ' Calcolo dei coefficienti necessari alla ' implementazione ricorsiva, multi frequenza ' dello spettro mobile su NsSpm campioni. ' Le frequenze, a cui calcolare gli spettri, ' devono essere definite nel vettore Freq(0 to NFreqSPM_1): ' Dim I&, OmegaT# ' NsSpm = NsSpm_I ' N° di campioni per il calcolo dello spettro mobile. NsSpm1 = NsSpm + 1 NFreqSPM_1 = NFRE - 1 ' N° di frequenze - 1 a cui calcolare gli spettri. NormSpm = CDbl(NsSpm) / 2# ' Fattore di normalizzazione sul N° di campioni. ' ReDim z1(0 To 0) ' z1(0) = e^j*omega*T. ReDim SmRE(0 To NsSpm) ' Registro a scorrimento dei campioni del segnale. ReDim SqRE(0 To NsSpm) ' Registro a scorrimento dei campioni del segnale al quadrato. VMed = 0# ' Valore medio corrente. ' ReDim Sn_MF_1(0 To NFreqSPM_1) ReDim z1_MF_1(0 To NFreqSPM_1) ReDim z1_MF_N(0 To NFreqSPM_1) ' For I = 0 To NFreqSPM_1 OmegaT = PI2 * Freq(I) / Fs ' Per frequenze in [Hz]. 'OmegaT = PI2 * Freq(I) ' Per frequenze in [f/fs]. z1(0) = CExp(CCmp(0#, OmegaT)) ' ' Per calcolo ricorsivo: Sn_MF_1(I) = CCmp(0#, 0#) z1_MF_1(I) = CPtN(z1(0), -1) z1_MF_N(I) = CPtN(z1(0), -NsSpm) Next I ' ' ' End Sub Public Sub SPM_MF_R(ByVal IMed As Long, ByVal Xn As Single, _ ByRef SpmMF() As Double, ByRef MedMF As Double) ' ' Routine per il calcolo ricorsivo dello spettro mobile (sliding ' spectrum) multi frequenza su NsSpm campioni. ' Devono essere stati definiti, precedentemente, i coefficienti ' z1_MF_1(), z1_MF_N() ed un vettore SmRE(0 To NsSpm) da usare come ' registro a scorrimento dei campioni del segnale. ' Ritorna anche il valore della media (o del valor efficace) mobile ' degli ultimi NsSpm campioni. ' In questa routine, per ragioni di velocita', le funzioni di numeri ' complessi sono state sostituite con il loro sviluppo in linea. ' Dim K& ' On Error Resume Next ' Gestisce gli errori di troncamento nel calcolo di MedMF. ' ' SmRE() usato come buffer circolare: Ic = ((Ic + 1) Mod NsSpm1) Kc = ((Ic + 1) Mod NsSpm1) SmRE(Ic) = CDbl(Xn) ' For K = 0 To NFreqSPM_1 'XZ1_N = CMol(CCmp(SmRE(Kc), 0#), z1_MF_N) Xz1_N.Re = SmRE(Kc) * z1_MF_N(K).Re Xz1_N.Im = SmRE(Kc) * z1_MF_N(K).Im 'SnI = CDif(CCmp(SmRE(Ic), 0#), XZ1_N) SnI.Re = SmRE(Ic) - Xz1_N.Re SnI.Im = -Xz1_N.Im ' 'YZ1_1 = CMol(Sn_MF_1(K), z1_MF_1(K)) Yz1_1.Re = Sn_MF_1(K).Re * z1_MF_1(K).Re - Sn_MF_1(K).Im * z1_MF_1(K).Im Yz1_1.Im = Sn_MF_1(K).Re * z1_MF_1(K).Im + Sn_MF_1(K).Im * z1_MF_1(K).Re 'Sn = CSom(SnI, YZ1_1) Sn.Re = SnI.Re + Yz1_1.Re Sn.Im = SnI.Im + Yz1_1.Im Sn_MF_1(K) = Sn ' 'SpmMF(K) = CAbs(Sn) / NormSpm SpmMF(K) = Sqr(Sn.Re * Sn.Re + Sn.Im * Sn.Im) / NormSpm Next K ' ' Calcolo la media/valor efficace mobile: If (NsSpm < IMed) Then ' Media mobile del valore assoluto del segnale: 'VMed = VMed + Abs(SmRE(Ic)) - Abs(SmRE(Kc)) 'MedMF = VMed / CDbl(NsSpm) ' ' Valore efficace mobile: SqRE(Ic) = SmRE(Ic) * SmRE(Ic) VMed = VMed + SqRE(Ic) - SqRE(Kc) MedMF = Sqr(VMed / CDbl(NsSpm)) ' Else ' Gestisco l' inizio del segnale: If (IMed = 1) Then VMed = 0# ' ' Media mobile del valore assoluto del segnale: 'VMed = VMed + Abs(SmRE(Ic)) 'MedMF = VMed / CDbl(IMed) ' ' Valore efficace mobile: SqRE(Ic) = SmRE(Ic) * SmRE(Ic) VMed = VMed + SqRE(Ic) MedMF = Sqr(VMed / CDbl(IMed)) End If ' ' ' End Sub Public Function SPM_R(ByVal Xn As Long) As Long 'Public Function SPM_R(ByVal Xn As Single) As Single ' ' Funzione per il calcolo ricorsivo dello spettro mobile ' (sliding spectrum) su NsSpm campioni. ' Devono essere stati definiti, precedentemente, i coefficienti ' Z1_1, Z1_N ed un vettore SmRE(0 To NsSpm) da usare come ' registro a scorrimento dei campioni del segnale. ' In questa routine, per ragioni di velocita', le funzioni di numeri ' complessi sono state sostituite con il loro sviluppo in linea. ' ' Versione per filtro del trigger di AudioCardDSP (24/04/2009). ' ' SmRE() usato come buffer circolare: Ic = ((Ic + 1) Mod NsSpm1) Kc = ((Ic + 1) Mod NsSpm1) SmRE(Ic) = Xn ' ' XZ1_N = CMol(CCmp(SmRE(NsSpm), 0#), Z1_N): Xz1_N.Re = SmRE(Kc) * z1_N.Re Xz1_N.Im = SmRE(Kc) * z1_N.Im ' ' SnI = CDif(CCmp(SmRE(0), 0#), XZ1_N): SnI.Re = SmRE(Ic) - Xz1_N.Re SnI.Im = -Xz1_N.Im ' ' Yz1_1 = CMol(Sn_1, z1_1): Yz1_1.Re = Sn_1.Re * z1_1.Re - Sn_1.Im * z1_1.Im Yz1_1.Im = Sn_1.Re * z1_1.Im + Sn_1.Im * z1_1.Re ' ' Sn = CSom(SnI, Yz1_1): Sn.Re = SnI.Re + Yz1_1.Re Sn.Im = SnI.Im + Yz1_1.Im Sn_1 = Sn ' ' SPM_R = CAbs(Sn) / NormSpm: SPM_R = Sqr(Sn.Re * Sn.Re + Sn.Im * Sn.Im) / NormSpm ' ' ' End Function Public Function MCorr_FT(D1() As Double, D2() As Double, _ ByVal R1 As Long, ByVal R2 As Long, _ Optional ByVal bUnbiased As Boolean = False, _ Optional ByVal bNoCC As Boolean = False) As Double() ' ' Espressione della mutua correlazione non circolare: ' MCorr_FT(R) = (L / Ld) * IFFT((1 / L) * Coniugata(FFT(D1(I))) * FFT(D2(I))) ' con: N1 = UBound(D1), N2 = UBound(D2), ' -N1 <= R1 <= 0, 0 <= R2 <= N2, R1 <= R <= R2, ' L = MIN(N1 + 1, N2 + 1), Ld = MIN(N1 + R, N2) - MAX(0, R) + 1 ' ' Ritorna un vettore contenente la funzione di mutua correlazione ' fra i dati in D1() e quelli in D2() calcolata su R1 --> R2 ritardi ' [i.e. MCorr_FT(R1 To R2)]. ' I segnali D1() e D2() sono considerati finiti. ' Se bUnbiased = True i valori della mutua correlazione vengono ' corretti con 1 / (N - r). ' Se bNoCC = True dai segnali viene eliminata la componente continua. ' ' Il calcolo viene effettuato con FFT e IFFT. ' Algoritmo tratto da: "Theory and Application of Digital ' Signal Processing" di L. Rabiner e B. Gold. - pg. 403 ' e corretto con le equazioni di "Digital Time Series Analysis" ' di R. K. Otnes, L. Enochson. - pg. 247, 248. ' ' Ver. 18/10/2005 modificata per AudioCardDSP. ' Dim M&, M2_1&, N1&, N2&, r& Dim LNorm#, dd1#(), dd2#(), NBytes& Dim Re1#(), Im1#(), Re2#(), Im2#(), s#(), f#() ' N1 = UBound(D1) N2 = UBound(D2) ' M = Ceil(Log((N1 + N2 + 1)) / Log2) M2_1 = (2 ^ M) - 1 ' ' I dati devono avere media zero? If bNoCC Then dd1() = EliminaCC(D1()) dd2() = EliminaCC(D2()) Else dd1() = D1() dd2() = D2() End If ' ReDim Preserve dd1(0 To M2_1) ' Aggiungo zeri per avere una potenza di 2. ReDim Preserve dd2(0 To M2_1) ' " " " " " " " " ' Sposto i valori per comporre ' opportunamente il secondo segnale: MoveMemory dd2(N1 + 1), dd2(0), 8 * (N2 + 1) ZeroMemory dd2(0), 8 * (N1 + 1) ' For R = N2 To 0 Step -1 ' dd2(N1 + R + 1) = dd2(R) ' dd2(R) = 0# ' Next R ' ' Calcolo la FFT dei due segnali: FFT_D2 dd1(), dd2(), Re1(), Im1(), Re2(), Im2(), M2_1 + 1 ' ' Moltiplico la seconda trasformata ' per la coniugata della prima: ReDim s(0 To M2_1), f(0 To M2_1) ' Uso S() e F() come vettori d' appoggio. For r = 0 To M2_1 s(r) = (Re1(r) * Re2(r) + Im1(r) * Im2(r)) ' Parte reale. f(r) = (Re1(r) * Im2(r) - Re2(r) * Im1(r)) ' Parte immaginaria. Next r ' ' Calcolo l' antitrasformata, corrispondente ' alla funzione di mutua correlazione: IFFT s(), f(), Re1(), Im1() ' Uso Re1() e Im1() come vettori d' appoggio. ' ' Sistemo i valori nell' ordine dei ritardi: ReDim Im1(R1 To R2) ' Uso Im1() come vettore d' appoggio. NBytes = 8 * (R2 - R1 + 1) ' I valori > N1 + N2 + 1 sono nulli e/o simmetrici. CopyMemory Im1(R1), Re1(N1 + R1 + 1), NBytes ' For R = R1 To R2 ' Im1(R) = Re1(N1 + R + 1) ' Next R ' If bUnbiased Then ' Correzione dei valori della mutua ' correlazione con 1 / (N - r): For r = R1 To R2 ' LNorm = MIN(N1 + R, N2) - MAX(0, R) + 1: If (N1 + r < N2) Then LNorm = CDbl(N1 + r + 1) Else LNorm = CDbl(N2 + 1) End If If (0 < r) Then LNorm = LNorm - CDbl(r) ' Im1(r) = Im1(r) / LNorm Next r ' Else ' Correzione dei valori della mutua ' correlazione con 1 / N: 'LNorm = CDbl(MIN0(N1 + 1, N2 + 1)) If (N1 <= N2) Then LNorm = CDbl(N1 + 1) Else LNorm = CDbl(N2 + 1) End If For r = R1 To R2 Im1(r) = Im1(r) / LNorm Next r End If ' MCorr_FT = Im1() ' ' ' End Function Public Function MAutoCorr_FT(D1() As Double, _ ByVal R1 As Long, ByVal R2 As Long, _ Optional ByVal bUnbiased As Boolean = True, _ Optional ByVal bNoCC As Boolean = True) As Double() ' ' Espressione dell' auto correlazione non circolare: ' MCorr_FT(R) = (L / Ld) * IFFT((1 / L) * Coniugata(FFT(D1(I))) * FFT(D1(I))) ' ' con: D1(0 To N1_C - 1), ' MM_C = Ceil(Log(CDbl(2 * N1_C - 1)) / Log2), ' N_C = 2 ^ MM_C, NFre_C = N_C / 2, ' -N1_C < R1 <= 0, 0 <= R2 < N1_C, R1 <= R <= R2, ' L = N1_C, Ld = N1_C - Abs(R) ' ' Ritorna un vettore contenente la funzione di auto correlazione ' dei dati in D1() calcolata su R1 --> R2 ritardi [i.e. MAutoCorr_FT(R1 To R2)]. ' ' I segnali D1() sono considerati finiti. ' Se bUnbiased = True i valori dell' auto correlazione vengono ' corretti con 1 / (N - |r|). ' Se bNoCC = True dai segnali viene eliminata la componente continua. ' ' Il calcolo viene effettuato con SFFTBF_Corr e SFFTBB_Corr e ' deve essere preceduto da una chiamata alla Sub SFFTBI_Corr. ' ' Algoritmo tratto da: "Theory and Application of Digital ' Signal Processing" di L. Rabiner e B. Gold. - pg. 403 ' e corretto con le equazioni di "Digital Time Series Analysis" ' di R. K. Otnes, L. Enochson. - pg. 247, 248. ' ' Ver. 18/10/2005 modificata per AudioCardDSP. ' Ver. 22/05/2007 con SFFTBF_Corr e SFFTBB_Corr. ' Dim r&, NBytes& Dim dd1#(), Sxx# Dim Corr#(), LNorm# ' ' I dati devono avere media zero? If bNoCC Then dd1() = EliminaCC(D1()) Else dd1() = D1() End If ' ReDim Preserve dd1(0 To N_C - 1) ' Aggiungo zeri per avere una potenza di 2. ' ' Calcolo la FFT del segnale: SFFTBF_Corr dd1(), Re1_C(), Im1_C() ' ' (A) - Calcolo la trasformata del segnale ' ritardato di N1_C campioni: 'For R = 0 To NFre_C ' Re2_C(R) = (Re1_C(R) * WnRe(R) - Im1_C(R) * WnIm(R)) ' Parte reale. ' Im2_C(R) = (WnRe(R) * Im1_C(R) + Re1_C(R) * WnIm(R)) ' Parte immaginaria. 'Next R ' ' (B) - Moltiplico la seconda trasformata ' per la coniugata della prima: 'For R = 0 To NFre_C ' S(R) = (Re1_C(R) * Re2_C(R) + Im1_C(R) * Im2_C(R)) ' Parte reale. ' F(R) = (Re1_C(R) * Im2_C(R) - Re2_C(R) * Im1_C(R)) ' Parte immaginaria. 'Next R ' ' Compatto i calcoli (A) e (B) precedenti per ' aumentarne la velocita' di esecuzione: For r = 0 To NFre_C Sxx = Re1_C(r) * Re1_C(r) + Im1_C(r) * Im1_C(r) s(r) = Sxx * WnRe(r) ' Parte reale. f(r) = Sxx * WnIm(r) ' Parte immaginaria. Next r ' ' Calcolo l' antitrasformata, corrispondente ' alla funzione di mutua correlazione: SFFTBB_Corr s(), f(), dd1() ' Uso dd1() come vettore d' appoggio. ' ' Sistemo i valori nell' ordine dei ritardi: ReDim Corr(R1 To R2) ' Uso Corr() come vettore d' appoggio. NBytes = 8 * (R2 - R1 + 1) ' I valori > 2 * N1_C - 1 sono nulli e/o simmetrici. CopyMemory Corr(R1), dd1(N1_C + R1), NBytes ' If bUnbiased Then ' Correzione dei valori della mutua ' correlazione con 1 / (N - |r|): For r = R1 To R2 LNorm = CDbl(N1_C - Abs(r)) Corr(r) = Corr(r) / LNorm Next r ' Else ' Correzione dei valori della mutua ' correlazione con 1 / N: LNorm = CDbl(N1_C) For r = R1 To R2 Corr(r) = Corr(r) / LNorm Next r End If ' MAutoCorr_FT = Corr() ' ' ' End Function Public Function MMutuaCorr_FT(D1() As Double, D2() As Double, _ ByVal R1 As Long, ByVal R2 As Long, _ Optional ByVal bUnbiased As Boolean = False, _ Optional ByVal bNoCC As Boolean = False) As Double() ' ' Espressione della mutua correlazione non circolare: ' MMutuaCorr_FT(R) = (L / Ld) * IFFT((1 / L) * Coniugata(FFT(D1(I))) * FFT(D2(I))) ' ' con: D1(0 To N1_C - 1), D2(0 To N2_C - 1), ' MM_C = Ceil(Log(CDbl(N1_C + N2_C - 1)) / Log2), ' N_C = 2 ^ MM_C, NFre_C = N_C / 2, ' -N1_C < R1 <= 0, 0 <= R2 < N2_C, R1 <= R <= R2, ' L = MIN(N1_C, N2_C), Ld = MIN(N1_C + R, N2_C) - MAX(0, R) ' ' Ritorna un vettore contenente la funzione di mutua correlazione ' fra i dati in D1() e quelli in D2() calcolata su R1 --> R2 ritardi ' [i.e. MAutoCorr_FT(R1 To R2)]. ' ' I segnali D1() e D2() sono considerati finiti. ' Se bUnbiased = True i valori dell' auto correlazione vengono ' corretti con 1 / (N - |r|). ' Se bNoCC = True dai segnali viene eliminata la componente continua. ' ' Il calcolo viene effettuato con SFFTBF_Corr e SFFTBB_Corr e ' deve essere preceduto da una chiamata alla Sub SFFTBI_Corr. ' ' Algoritmo tratto da: "Theory and Application of Digital ' Signal Processing" di L. Rabiner e B. Gold. - pg. 403 ' e corretto con le equazioni di "Digital Time Series Analysis" ' di R. K. Otnes, L. Enochson. - pg. 247, 248. ' ' Ver. 18/10/2005 modificata per AudioCardDSP. ' Ver. 22/05/2007 con SFFTBF_Corr e SFFTBB_Corr. ' ' Dim r&, NBytes& Dim dd1#(), dd2#() Dim Corr#(), LNorm# ' ' I dati devono avere media zero? If bNoCC Then dd1() = EliminaCC(D1()) dd2() = EliminaCC(D2()) Else dd1() = D1() dd2() = D2() End If ' ReDim Preserve dd1(0 To N_C - 1) ' Aggiungo zeri per avere una potenza di 2. ReDim Preserve dd2(0 To N_C - 1) ' Aggiungo zeri per avere una potenza di 2. ' Sposto i valori per comporre ' opportunamente il secondo segnale ' ritardato di N1_C campioni: MoveMemory dd2(N1_C), dd2(0), 8 * N2_C ZeroMemory dd2(0), 8 * N1_C ' ' Calcolo le FFT dei segnali: SFFTBF_Corr dd1(), Re1_C(), Im1_C() SFFTBF_Corr dd2(), Re2_C(), Im2_C() ' ' Moltiplico la seconda trasformata ' per la coniugata della prima: For r = 0 To NFre_C s(r) = (Re1_C(r) * Re2_C(r) + Im1_C(r) * Im2_C(r)) ' Parte reale. f(r) = (Re1_C(r) * Im2_C(r) - Re2_C(r) * Im1_C(r)) ' Parte immaginaria. Next r ' ' Calcolo l' antitrasformata, corrispondente ' alla funzione di mutua correlazione: SFFTBB_Corr s(), f(), dd1() ' Uso dd1() come vettore d' appoggio. ' ' Sistemo i valori nell' ordine dei ritardi: ReDim Corr(R1 To R2) ' Uso Corr() come vettore d' appoggio. NBytes = 8 * (R2 - R1 + 1) ' I valori > N1_C + N2_C - 1 sono nulli e/o simmetrici. CopyMemory Corr(R1), dd1(N1_C + R1), NBytes ' If bUnbiased Then ' Correzione dei valori della mutua ' correlazione con 1 / (N - |r|): For r = R1 To R2 ' LNorm = MIN0(N1_C + R, N2_C) - MAX0(0, R): If (N1_C + r < N2_C) Then LNorm = CDbl(N1_C + r) Else LNorm = CDbl(N2_C) End If If (0 < r) Then LNorm = LNorm - CDbl(r) ' Corr(r) = Corr(r) / LNorm Next r ' Else ' Correzione dei valori della mutua ' correlazione con 1 / N: ' LNorm = CDbl(MIN0(N1_C, N2_C)): If (N1_C <= N2_C) Then LNorm = CDbl(N1_C) Else LNorm = CDbl(N2_C) End If ' For r = R1 To R2 Corr(r) = Corr(r) / LNorm Next r End If ' MMutuaCorr_FT = Corr() ' ' ' End Function Public Sub FFT_D2(D1() As Double, D2() As Double, _ R1() As Double, X1() As Double, R2() As Double, X2() As Double, ByRef NVAL As Long) ' ' Definizioni: ' ' FFT = Fast Fourier Transform. ' ' Entra con: ' D1() = Vettore dei dati Reali del primo Segnale(t). ' D2() = Vettore dei dati Reali del secondo Segnale(t). ' NVAL = Numero max. dei dati in D1(), D2() da usare. ' ' Esce con: ' R1() = Vettore dei valori Reali della Trasformata(f) del primo segnale. ' X1() = Vettore dei valori Immaginari della Trasformata(f) del primo segnale. ' R2() = Vettore dei valori Reali della Trasformata(f) del secondo segnale. ' X2() = Vettore dei valori Immaginari della Trasformata(f) del secondo segnale. ' NVAL = Numero dei dati usati dalla Trasformata(f). ' NFRE = Numero dei valori di Frequenza calcolati. ' ' Calcola, contemporaneamente, le trasformate di due segnali ' reali definiti in D1(0 To NVAL - 1) e D2(0 To NVAL - 1). ' ' I valori delle Trasformate in R1(), X1(), R2(), X2() NON vengono ' normalizzati sulla lunghezza del Segnale 2^M. ' ' I vettori R1(), X1() e R2(), X2() sono dimensionati in questa routine. ' Gli NN dati nei vettori sono organizzati come Vettore(0 To NN - 1). ' ' Tradotta in Basic da: "Theory and Application of Digital ' "Theory and Application of Digital Signal Processing" ' di L. Rabiner e B. Gold. - pg. 367. ' e con le formule di: ' "Digital Time Series Analysis" ' di R. Otnes e L. Enochson - pg. 175. ' ' Ver. 18/10/2005 modificata per AudioCardDSP. ' Dim M&, I&, NMN1&, NFRE&, J&, K&, L&, LE&, LE1&, IP&, JF&, NBytes& Dim TSwap#, Ur#, Ui#, wr#, wi#, Tr#, Ti# ' M = Int(Log(CDbl(NVAL) + 0.5) / Log2) NVAL = 2 ^ M NFRE = NVAL / 2 NMN1 = NVAL - 1 ' ReDim Re#(0 To NVAL - 1), Im#(0 To NVAL - 1) ReDim R1(0 To NVAL - 1), R2(0 To NVAL - 1), X1(0 To NVAL - 1), X2(0 To NVAL - 1) ' ' Compongo i vettori Re(), Im() da trasformare ' con i dati reali in ingresso: NBytes = 8 * NVAL CopyMemory Re(0), D1(0), NBytes CopyMemory Im(0), D2(0), NBytes ' For I = 0 To NVAL - 1 ' Re(I) = D1(I) ' Im(I) = D2(I) ' Next I ' J = 0 For I = 0 To NMN1 - 1 If (I < J) Then ' DSWAP Re(J), Re(I): TSwap = Re(J) Re(J) = Re(I) Re(I) = TSwap ' DSWAP Im(J), Im(I) TSwap = Im(J) Im(J) = Im(I) Im(I) = TSwap End If K = NFRE Do While (K - 1 < J) J = J - K K = K / 2 Loop J = J + K Next I ' For L = 1 To M LE = 2 ^ L LE1 = LE / 2 Ur = 1# Ui = 0# wr = Cos(PI / LE1) wi = Sin(PI / LE1) For J = 0 To LE1 - 1 For I = J To NVAL - 1 Step LE IP = I + LE1 Tr = Re(IP) * Ur - Im(IP) * Ui Ti = Re(IP) * Ui + Im(IP) * Ur Re(IP) = Re(I) - Tr Im(IP) = Im(I) - Ti Re(I) = Re(I) + Tr Im(I) = Im(I) + Ti Next I Tr = Ur Ti = Ui Ur = Tr * wr - Ti * wi Ui = Tr * wi + Ti * wr Next J Next L ' ' Calcolo delle parti reali ed immaginarie delle trasformate: R1(0) = Re(0) X1(0) = 0# R2(0) = Im(0) X2(0) = 0# ' For JF = 1 To NVAL - 1 R1(JF) = (Re(JF) + Re(NVAL - JF)) / 2# X1(JF) = (Im(JF) - Im(NVAL - JF)) / 2# R2(JF) = (Im(JF) + Im(NVAL - JF)) / 2# X2(JF) = (Re(NVAL - JF) - Re(JF)) / 2# Next JF ' ' ' End Sub Public Sub IFFT(r() As Double, X() As Double, Dr() As Double, dx() As Double, _ Optional ByRef NVAL As Long) ' ' Definizioni: ' ' IFFT = Inverse Fast Fourier Transform. ' ' Entra con: ' R() = Vettore dei valori Reali della Trasformata(f). ' X() = Vettore dei valori Immaginari della Trasformata(f). ' ' Esce con: ' Dr() = Vettore dei dati Reali del Segnale(t). ' Dx() = Vettore dei dati Immaginari del Segnale(t). ' NVAL = Numero dei dati usati per l' antitrasformata. ' ' I vettori Dr(), Dx() sono dimensionati in questa routine. ' Gli NN dati nei vettori sono organizzati come Vettore(0 To NN - 1). ' ' Algoritmo tratto da: "Theory and Application of Digital ' Signal Processing" di L. Rabiner e B. Gold. - pg. 371. ' Dim NFRE&, M&, I&, NMN1&, J& Dim K&, L&, LE&, LE1&, IP&, JF& Dim Ur#, Ui#, wr#, wi#, Tr#, Ti# ' NVAL = UBound(r) + 1 M = Int(Log(CDbl(NVAL) + 0.5) / Log2) NVAL = 2 ^ M NFRE = NVAL / 2 NMN1 = NVAL - 1 ' ReDim Dr(0 To NVAL - 1), dx(0 To NVAL - 1), d(0 To NVAL - 1) ' ' Se non si vogliono sfruttare le proprieta' ' "transform-in-place" di questa implementazione: For I = 0 To NVAL - 1 Dr(I) = r(I) / CDbl(NVAL) dx(I) = -X(I) / CDbl(NVAL) Next I ' J = 0 For I = 0 To NMN1 - 1 If (I < J) Then DSWAP Dr(J), Dr(I) DSWAP dx(J), dx(I) End If K = NFRE Do While (K - 1 < J) J = J - K K = K / 2 Loop J = J + K Next I ' For L = 1 To M LE = 2 ^ L LE1 = LE / 2 Ur = 1# Ui = 0# wr = Cos(PI / LE1) wi = Sin(PI / LE1) For J = 0 To LE1 - 1 For I = J To NVAL - 1 Step LE IP = I + LE1 Tr = Dr(IP) * Ur - dx(IP) * Ui Ti = Dr(IP) * Ui + dx(IP) * Ur Dr(IP) = Dr(I) - Tr dx(IP) = dx(I) - Ti Dr(I) = Dr(I) + Tr dx(I) = dx(I) + Ti Next I Tr = Ur Ti = Ui Ur = Tr * wr - Ti * wi Ui = Tr * wi + Ti * wr Next J Next L ' ' ' End Sub Public Function MCorr(D1() As Double, D2() As Double, _ ByVal R1 As Long, ByVal R2 As Long, _ Optional ByVal bUnbiased As Boolean = False, _ Optional ByVal bNoCC As Boolean = False) As Double() ' ' Espressione della mutua correlazione non circolare: ' I2 ' MCorr(R) = (1 / (I2 - I1 + 1)) * Som(D1(I) * D2(I + R)) con: N1 = UBound(D1) ' I=I1 " N2 = UBound(D2) ' " -N1 <= R1 < R2 ' " R1 < R2 <= N2 ' " R1 <= R <= R2 ' " I1 = MAX(0, -R) ' " I2 = MIN(N1, N2 - R) ' ' Ritorna un vettore contenente la funzione di mutua correlazione ' fra i dati in D1() e quelli in D2() calcolata su R1 --> R2 ritardi ' [i.e. MCorr(R1 To R2)]. ' Il segnale D2() e' considerato finito e costituito da UBound(D2) + 1 ' campioni. ' Se bUnbiased = True i valori della mutua correlazione vengono ' corretti con 1 / (N - r). ' Se bNoCC = True dai segnali viene eliminata la componente continua. ' ' Il calcolo viene effettuato per convoluzione diretta ' e corretto con le equazioni di "Digital Time Series Analysis" ' di R. K. Otnes, L. Enochson. - pg. 247. ' ' Ver. 16/10/2005 modificata per AudioCardDSP. ' Dim I&, I1&, I2&, N1&, N2&, r& ', R1&, R2& Dim LNorm#, dd1#(), dd2#() ' N1 = UBound(D1) N2 = UBound(D2) ' ReDim RR#(R1 To R2) ReDim MCorr(R1 To R2) LNorm = CDbl(MIN0(N1 + 1, N2 + 1)) ' ' I dati devono avere media zero? If bNoCC Then dd1() = EliminaCC(D1()) dd2() = EliminaCC(D2()) Else dd1() = D1() dd2() = D2() End If ' ' Calcolo la funzione di mutua correlazione: For r = R1 To R2 I1 = MAX0(0, -r) I2 = MIN0(N1, N2 - r) For I = I1 To I2 RR(r) = RR(r) + dd1(I) * dd2(I + r) Next I ' If bUnbiased Then RR(r) = RR(r) / CDbl(I2 - I1 + 1) Else RR(r) = RR(r) / LNorm End If Next r ' MCorr = RR() ' ' ' End Function Public Function IIR(ByVal VIn As Single) As Single ' ' Funzione di filtrazione IIR dei segnali. ' Devono essere state definite, precedentemente, le matrici ' Ac(0 To NK, 1 To NCel) e Bc(0 To NK, 1 To NCel) contenenti ' i coefficienti del filtro ed una matrice W(0 To NK, 1 To NCel) ' da usare come registro a scorrimento dei campioni da filtrare. ' ' Codice tradotto ed adattato dal FORTRAN di: ' "Digital Filters and their Applications" ' di V. Cappellini, A. G. Constantinides, P. Emiliani. ' pg. 373. ' ' Routine modificata per strumenti di misura. ' ------------------------------------------- ' Dim C&, K&, Y! ' Y = VIn For C = 1 To NCel w(0, C) = Y Y = 0! For K = NK To 1 Step -1 w(0, C) = w(0, C) - Bc(K, C) * w(K, C) Y = Y + Ac(K, C) * w(K, C) w(K, C) = w(K - 1, C) Next K Y = Y + Ac(0, C) * w(0, C) Next C ' IIR = Y ' ' ' End Function Public Function FIR(ByVal VIn As Single) As Single ' ' Funzione di filtrazione FIR dei segnali. ' Devono essere stati definiti, precedentemente, un vettore ' Hc(0 To NCoeff - 1) contenente gli NCoeff coefficienti del ' filtro ed un vettore SRE(-1 To NCoeff - 1) da usare come ' registro a scorrimento dei campioni da filtrare. ' ' Codice tradotto ed adattato dal FORTRAN di: ' "Digital Filters and their Applications" ' di V. Cappellini, A. G. Constantinides, P. Emiliani. ' pg. 371. ' ' Routine modificata per strumenti di misura. ' ------------------------------------------- ' Dim K&, VOut! ' SRE(-1) = VIn For K = UBound(SRE) To 0 Step -1 SRE(K) = SRE(K - 1) VOut = VOut + SRE(K) * Hc(K) Next K ' FIR = VOut ' ' ' End Function Public Sub SFFTBB(ByVal N As Long, ByVal NFRE As Long, _ Re() As Double, Im() As Double, ByRef d() As Double) ' ' SFFTBB( X, N, MM, S1, C1, S3, C3, ITAB ) ' ' A real-valued, in place, split-radix IFFT program ' Hermitian symmetric input and real output in array X ' Length is N = 2 ** MM ' Decimation-in-frequency, cos/sin in second loop with table look-up ' Input order: ' [ Re(0), Re(1), ..., Re(N/2), Im(N/2-1), ..., Im(1) ] ' ' S1 - array of sin() table (length >= N/8-1) ' C1 - array of cos() table (length >= N/8-1) ' S3 - array of sin() table (length >= N/8-1) ' C3 - array of cos() table (length >= N/8-1) ' ITAB - integer bit reversal table (length >= sqrt(2n) ) ' ' The initialization routine SFFTBI must be called prior to calling ' this routine. SFFTBI need not be called again unless N changes. ' ' Original code (IRVFFT) written by H.V. Sorensen, ' Rice University, Oct. 1985 ' ' Modifications made by Steve Kifowit, 26 June 1997 ' -- table look-up of sines and cosines ' -- incorporation of bit reversal table ' -- quick return ' ' Tradotta dal FORTRAN e modificata da F. Languasco 15/10/2005. ' ' Entra con: ' N = Numero max. di dati in D() da calcolare. ' NFRE = Numero dei valori di Frequenza nei vettori Re() e Im(). ' Re() = Vettore dei valori Reali della Trasformata(f). ' Im() = Vettore dei valori Immaginari della Trasformata(f). ' ' Esce con: ' D() = Vettore dei valori Reali del Segnale(t), ridimensionato ' in questa routine. ' ' Gli NFRE + 1 dati nei vettori Re() e Im() sono organizzati come ' Vettore(0 To NFRE); gli N dati nel vettore D() sono organizzati ' come D(0 To N - 1). ' Dim J&, I&, K&, Ic&, ID&, I0&, I1&, I2&, I3&, I4&, I5&, I6&, I7&, I8& Dim N1&, N2&, N4&, N8&, NN&, It& Dim XT#, R1#, t1#, T2#, T3#, T4#, T5# Dim CC1#, SS1#, CC3#, SS3# ReDim d#(0 To N - 1) ' If (N = 1) Then Exit Sub ' ' Sistemo i valori dei vettori Re(), Im() in D1(): ' For I = 0 To NFre ' D1(I + 1) = Re(I) ' Next I CopyMemory D1(1), Re(0), 8 * (NFRE + 1) For I = 1 To NFRE - 1 D1(N - I + 1) = Im(I) Next I ' N2 = 2 * N NN = 1 For K = 1 To MM - 1 Ic = 0 ID = N2 N2 = N2 / 2 N4 = N2 / 4 N8 = N4 / 2 17 For I = Ic To N - 1 Step ID I1 = I + 1 I2 = I1 + N4 I3 = I2 + N4 I4 = I3 + N4 t1 = D1(I1) - D1(I3) D1(I1) = D1(I1) + D1(I3) D1(I2) = 2 * D1(I2) D1(I3) = t1 - 2 * D1(I4) D1(I4) = t1 + 2 * D1(I4) If (N4 = 1) Then GoTo 15 I1 = I1 + N8 I2 = I2 + N8 I3 = I3 + N8 I4 = I4 + N8 t1 = (D1(I2) - D1(I1)) / Sqr2 T2 = (D1(I4) + D1(I3)) / Sqr2 D1(I1) = D1(I1) + D1(I2) D1(I2) = D1(I4) - D1(I3) D1(I3) = 2 * (-T2 - t1) D1(I4) = 2 * (-T2 + t1) 15 ' CONTINUE Next I Ic = 2 * ID - N2 ID = 4 * ID If (Ic < N - 1) Then GoTo 17 For J = 2 To N8 It = (J - 1) * NN CC1 = C1(It) SS1 = S1(It) CC3 = C3(It) SS3 = S3(It) Ic = 0 ID = 2 * N2 40 For I = Ic To N - 1 Step ID I1 = I + J I2 = I1 + N4 I3 = I2 + N4 I4 = I3 + N4 I5 = I + N4 - J + 2 I6 = I5 + N4 I7 = I6 + N4 I8 = I7 + N4 t1 = D1(I1) - D1(I6) D1(I1) = D1(I1) + D1(I6) T2 = D1(I5) - D1(I2) D1(I5) = D1(I2) + D1(I5) T3 = D1(I8) + D1(I3) D1(I6) = D1(I8) - D1(I3) T4 = D1(I4) + D1(I7) D1(I2) = D1(I4) - D1(I7) T5 = t1 - T4 t1 = t1 + T4 T4 = T2 - T3 T2 = T2 + T3 D1(I3) = T5 * CC1 + T4 * SS1 D1(I7) = -T4 * CC1 + T5 * SS1 D1(I4) = t1 * CC3 - T2 * SS3 D1(I8) = T2 * CC3 + t1 * SS3 30 ' CONTINUE Next I Ic = 2 * ID - N2 ID = 4 * ID If (Ic < N - 1) Then GoTo 40 20 ' CONTINUE Next J NN = 2 * NN Next K ' Ic = 1 ID = 4 70 For I0 = Ic To N Step ID I1 = I0 + 1 R1 = D1(I0) D1(I0) = R1 + D1(I1) D1(I1) = R1 - D1(I1) 60 ' CONTINUE Next I0 Ic = 2 * ID - 1 ID = 4 * ID If (Ic < N) Then GoTo 70 ' N1 = ITAB(1) For K = 2 To N1 I0 = N1 * ITAB(K) + 1 I = K J = I0 For It = 2 To ITAB(K) + 1 XT = D1(I) D1(I) = D1(J) D1(J) = XT I = I + N1 J = I0 + ITAB(It) 101 ' CONTINUE Next It 100 ' CONTINUE Next K ' For I = 1 To N D1(I) = D1(I) / N 99 ' CONTINUE Next I ' ' Sposto i dati calcolati nel vettore D(): CopyMemory d(0), D1(1), 8 * N ' ' ... End of subroutine SFFTBB ... ' End Sub Public Sub SFFTBB_Corr(Re_C() As Double, Im_C() As Double, ByRef d() As Double) ' ' SFFTBB( X, N_C, MM_C, S1_C, C1_C, S3_C, C3_C, ITAB_C ) ' ' A real-valued, in place, split-radix IFFT program ' Hermitian symmetric input and real output in array X ' Length is N_C = 2 ** MM_C ' Decimation-in-frequency, cos/sin in second loop with table look-up ' Input order: ' [ Re_C(0), Re_C(1), ..., Re_C(N_C/2), Im_C(N_C/2-1), ..., Im_C(1) ] ' ' S1_C - array of sin() table (length >= N_C/8-1) ' C1_C - array of cos() table (length >= N_C/8-1) ' S3_C - array of sin() table (length >= N_C/8-1) ' C3_C - array of cos() table (length >= N_C/8-1) ' ITAB_C - integer bit reversal table (length >= sqrt(2n) ) ' ' The initialization routine SFFTBI_Corr must be called prior to calling ' this routine. SFFTBI_Corr need not be called again unless N1_C or ' N2_C change. ' ' Original code (IRVFFT) written by H.V. Sorensen, ' Rice University, Oct. 1985 ' ' Modifications made by Steve Kifowit, 26 June 1997 ' -- table look-up of sines and cosines ' -- incorporation of bit reversal table ' -- quick return ' ' Tradotta dal FORTRAN e modificata da F. Languasco 15/10/2005. ' ' Entra con: ' N_C = Numero max. di dati in D() da calcolare. ' NFre_C = Numero dei valori di Frequenza nei vettori Re_C() e Im_C(). ' Re_C() = Vettore dei valori Reali della Trasformata(f). ' Im_C() = Vettore dei valori Immaginari della Trasformata(f). ' ' Esce con: ' D() = Vettore dei valori Reali del Segnale(t), ridimensionato ' in questa routine. ' ' Gli NFre_C + 1 dati nei vettori Re_C() e Im_C() sono organizzati come ' Vettore(0 To NFre_C); gli N_C dati nel vettore D() sono organizzati ' come D(0 To N_C - 1). ' ' Ver: 22/05/2007 per MAutoCorr_FT e MMutuaCorr_FT. ' Dim J&, I&, K&, Ic&, ID&, I0&, I1&, I2&, I3&, I4&, I5&, I6&, I7&, I8& Dim N1&, N2&, N4&, N8&, NN&, It& Dim XT#, R1#, t1#, T2#, T3#, T4#, T5# Dim CC1#, SS1#, CC3#, SS3# ReDim d#(0 To N_C - 1) ' If (N_C = 1) Then Exit Sub ' ' Sistemo i valori dei vettori Re_C(), Im_C() in D1_C(): CopyMemory D1_C(1), Re_C(0), 8 * (NFre_C + 1) For I = 1 To NFre_C - 1 D1_C(N_C - I + 1) = Im_C(I) Next I ' N2 = 2 * N_C NN = 1 For K = 1 To MM_C - 1 Ic = 0 ID = N2 N2 = N2 / 2 N4 = N2 / 4 N8 = N4 / 2 17 For I = Ic To N_C - 1 Step ID I1 = I + 1 I2 = I1 + N4 I3 = I2 + N4 I4 = I3 + N4 t1 = D1_C(I1) - D1_C(I3) D1_C(I1) = D1_C(I1) + D1_C(I3) D1_C(I2) = 2 * D1_C(I2) D1_C(I3) = t1 - 2 * D1_C(I4) D1_C(I4) = t1 + 2 * D1_C(I4) If (N4 = 1) Then GoTo 15 I1 = I1 + N8 I2 = I2 + N8 I3 = I3 + N8 I4 = I4 + N8 t1 = (D1_C(I2) - D1_C(I1)) / Sqr2 T2 = (D1_C(I4) + D1_C(I3)) / Sqr2 D1_C(I1) = D1_C(I1) + D1_C(I2) D1_C(I2) = D1_C(I4) - D1_C(I3) D1_C(I3) = 2 * (-T2 - t1) D1_C(I4) = 2 * (-T2 + t1) 15 ' CONTINUE Next I ' Ic = 2 * ID - N2 ID = 4 * ID If (Ic < N_C - 1) Then GoTo 17 ' For J = 2 To N8 It = (J - 1) * NN CC1 = C1_C(It) SS1 = S1_C(It) CC3 = C3_C(It) SS3 = S3_C(It) Ic = 0 ID = 2 * N2 40 For I = Ic To N_C - 1 Step ID I1 = I + J I2 = I1 + N4 I3 = I2 + N4 I4 = I3 + N4 I5 = I + N4 - J + 2 I6 = I5 + N4 I7 = I6 + N4 I8 = I7 + N4 t1 = D1_C(I1) - D1_C(I6) D1_C(I1) = D1_C(I1) + D1_C(I6) T2 = D1_C(I5) - D1_C(I2) D1_C(I5) = D1_C(I2) + D1_C(I5) T3 = D1_C(I8) + D1_C(I3) D1_C(I6) = D1_C(I8) - D1_C(I3) T4 = D1_C(I4) + D1_C(I7) D1_C(I2) = D1_C(I4) - D1_C(I7) T5 = t1 - T4 t1 = t1 + T4 T4 = T2 - T3 T2 = T2 + T3 D1_C(I3) = T5 * CC1 + T4 * SS1 D1_C(I7) = -T4 * CC1 + T5 * SS1 D1_C(I4) = t1 * CC3 - T2 * SS3 D1_C(I8) = T2 * CC3 + t1 * SS3 30 ' CONTINUE Next I ' Ic = 2 * ID - N2 ID = 4 * ID If (Ic < N_C - 1) Then GoTo 40 20 ' CONTINUE Next J NN = 2 * NN Next K ' Ic = 1 ID = 4 70 For I0 = Ic To N_C Step ID I1 = I0 + 1 R1 = D1_C(I0) D1_C(I0) = R1 + D1_C(I1) D1_C(I1) = R1 - D1_C(I1) 60 ' CONTINUE Next I0 ' Ic = 2 * ID - 1 ID = 4 * ID If (Ic < N_C) Then GoTo 70 ' N1 = ITAB_C(1) For K = 2 To N1 I0 = N1 * ITAB_C(K) + 1 I = K J = I0 For It = 2 To ITAB_C(K) + 1 XT = D1_C(I) D1_C(I) = D1_C(J) D1_C(J) = XT I = I + N1 J = I0 + ITAB_C(It) 101 ' CONTINUE Next It 100 ' CONTINUE Next K ' For I = 1 To N_C D1_C(I) = D1_C(I) / N_C 99 ' CONTINUE Next I ' ' Sposto i dati calcolati nel vettore D(): CopyMemory d(0), D1_C(1), 8 * N_C ' ' ... End of subroutine SFFTBB ... ' End Sub Public Function SintesiFIR_WM(ByVal Tipo$, ByVal NFilt As Long, _ ByVal F1 As Double, ByVal F2 As Double, ByVal WNome$, ByVal Par As Double) As Boolean ' ' Sintesi del filtro con Window Method. ' ' Routine modificata per strumenti di misura. ' ------------------------------------------- ' ' Parametri in ingresso. ' Tipo$: Tipo di filtro: ' "Low Pass" ' "High Pass" ' "Band Pass" ' "Band Stop" ' "Differentiator" ' "Hilbert Trans." ' NFilt: Numero di Coefficienti del filtro. ' F1: Prima Frequenza di Taglio (0 <= F1 < Fs/2). ' F2: Seconda Frequenza di Taglio (F1 < F2 <= Fs/2). ' WNome$: Nome del tipo di Window desiderato (vedi Function WinTipi()). ' Par: Parametro per certi tipi di Windows. ' ' I coefficienti del filtro vengono calcolati ' nel vettore Hc(0 To NFilt - 1). ' Dim I&, K&, KK&, NOdd&, NFilt1& Dim AA#, Ak#, Af1#, Af2#, Den# ReDim Hc(0 To NFilt - 1), SRE(-1 To NFilt - 1) ' On Error GoTo SintesiFIR_WM_ERR ' NOdd = NFilt Mod 2 NFilt1 = Int(NFilt / 2) AA = 0.5 * (1# - CDbl(NOdd)) ' WF() = WinProf(WNome$, NFilt, Par) ' ' Impulse response evaluation ' --------------------------- ' If NOdd = 1 Then Select Case Tipo$ Case "Low Pass" Hc(NFilt1) = 2# * F1 Case "High Pass" Hc(NFilt1) = 1# - 2# * F1 Case "Band Pass" Hc(NFilt1) = 2# * (F2 - F1) Case "Band Stop" Hc(NFilt1) = 1# - 2# * (F2 - F1) Case "Differentiator" Hc(NFilt1) = 0# Case "Hilbert Trans." Hc(NFilt1) = 0# End Select End If ' For I = 1 To NFilt1 K = NFilt1 - I Ak = CDbl(I) - AA Af1 = PI2 * Ak * F1 Af2 = PI2 * Ak * F2 Den = PI * Ak ' Select Case Tipo$ Case "Low Pass" Hc(K) = Sin(Af1) / Den ' Case "High Pass" Hc(K) = (Sin(PI * Ak) - Sin(Af1)) / Den ' Case "Band Pass" Hc(K) = (Sin(Af2) - Sin(Af1)) / Den ' Case "Band Stop" Hc(K) = (Sin(PI * Ak) - Sin(Af2) + Sin(Af1)) / Den ' Case "Differentiator" Hc(K) = -(Sin(Af1) - Af1 * Cos(Af1)) / (Den * Den) ' Case "Hilbert Trans." Hc(K) = -(Cos(Af1) - Cos(Af2)) / Den End Select ' ' Windowing: Hc(K) = Hc(K) * WF(K) Next I ' ' For K = 0 To NFilt1 - 1 KK = NFilt - 1 - K Hc(KK) = Switch(Tipo$ = "Differentiator", -Hc(K), _ Tipo$ = "Hilbert Trans.", -Hc(K), _ True, Hc(K)) Next K ' ' SintesiFIR_WM_ERR: SintesiFIR_WM = (Err.Number = 0) If (Err.Number <> 0) Then Dim M$ M$ = FiltroErrFIR$ & vbNewLine & vbNewLine M$ = M$ & "Errore " & Str$(Err.Number) & vbNewLine M$ = M$ & Err.Description MsgBox M$, vbCritical, " modDSP: SintesiFIR_WM" End If ' ' ' End Function Public Function SintesiIIR_BT(ByVal Trans$, ByVal Proty$, ByRef FCt As Double, _ Optional ByVal FAt As Double, Optional ByVal Att As Double, Optional ByVal Rip As Double, _ Optional ByVal F1 As Double, Optional ByVal F2 As Double, _ Optional ByVal q As Double, Optional ByRef AttTot As Double, _ Optional ByRef Alfa As Double) As Boolean ' ' Sintesi del filtro richiesto con il metodo ' di Trasformazione Bilineare. ' ' Routine modificata per strumenti di misura. ' ------------------------------------------- ' ' Parametri in ingresso. ' Trans$: Tipo di Trasformazione richiesta: ' "No Transformation" JBT = 0 ' "Low Pass -> Low Pass" JBT = 1 ' "Low Pass -> High Pass" JBT = 2 ' "Low Pass -> Band Pass" JBT = 3 ' "Low Pass -> Band Stop" JBT = 4 ' Filtri elementari: ' "RC Type Low Pass Filter" JBT = 5 ' "RC Type High Pass Filter" JBT = 6 ' "Resonator" JBT = 7 ' "Notch Filter" JBT = 8 ' "U(t) = a·I(t) + (1-a)·U(t-1)" JBT = 9 ' Proty$: Tipo di prototipo Passa Basso: ' "Butterworth" ' "Chebyshev" ' FCt: Frequenza di taglio (con atten. di -3 dB) ' o di risonanza (solo per JBT = 7 e JBT = 8). ' FAt: Per JBT < 5 e' la frequenza alla quale ' l' attenuazione deve essere di almeno Att# dB. ' Att: Attenuazione desiderata in [dB] alla frequenza FAt#. ' Rip: Ondulazione in Banda Passante espressa come ' percentuale del guadagno (solo per Chebyshev). ' F1: Prima frequenza di taglio dopo la trasformazione. ' F2: Seconda frequenza di taglio dopo la trasformazione. ' Q: Solo per JBT = 7 e JBT = 8, e' il Q del risonatore. ' Alfa: Alfa dato per filtro JBT = 9; per specificare, invece, ' la FCt passare Alfa = -1. ' ' Parametri in uscita ' AttTot: Attenuazione ottenuta in [dB] alla frequenza FAt#. ' FCt: Frequenza di taglio ottenuta per filtro JBT = 9 con Alfa dato. ' Alfa: Alfa ottenuto per filtro JBT = 9 con FCt dato. ' ' I coefficienti del filtro vengono calcolati nelle ' matrici Ac(0 To NK, 1 To NCel) e Bc(0 To NK, 1 To NCel) ' con NK = ordine del filtro e NCel = numero di sezioni. ' In questa versione viene anche dimensionata la matrice ' W() dei registri delle sezioni del filtro. ' Dim N&, I&, K&, NC&, M$ Dim Om#, Cm#, a#, A0#, A1#, A2#, B0#, B1#, B2# ' Variabili per filtri elementari. Dim f#, FF#, FCTA#, FATA#, Eps#, ATL#, res# ' Variabili per Dim AB#, AT#, Bt#, TETAP#, P1#, P2#, P3# ' filtri di tipo Dim CAPPA#, COSTA#, TOT#, CC#, AL#, Rn# ' Butterworth e Chebyshev. ' On Error GoTo SintesiIIR_BT_ERR ' '-- Filtri elementari ------------------------------------------------------------------ If (Trans$ = "RC Type Low Pass Filter") _ Or (Trans$ = "RC Type High Pass Filter") _ Or (Trans$ = "U(t) = a·I(t) + (1-a)·U(t-1)") Then ' Calcolo di tipi di filtro elementari del 1° ordine. ' Vedi Sub zNote: Om = PI2 * FCt ' Pulsazione di taglio [rad/s]. ' Select Case Trans$ Case "RC Type Low Pass Filter" Cm = Om / Tan(PI * FCt) ' Costante di mappatura. a = 1# + Cm / Om A0 = 1# / a A1 = 1# / a B0 = 1# B1 = (1# - Cm / Om) / a ' Case "RC Type High Pass Filter" Cm = Om / Tan(PI * FCt) ' Costante di mappatura. a = 1# + Cm / Om A0 = (Cm / Om) / a A1 = -(Cm / Om) / a B0 = 1# B1 = (1# - Cm / Om) / a ' Case "U(t) = a·I(t) + (1-a)·U(t-1)" ' a := Alfa If (Alfa < 0#) Then ' E' data FCt; calcola Alfa per avere la FCt voluta: Alfa = Cos(Om) - 1# + Sqr(Cos(Om) ^ 2 - 4# * Cos(Om) + 3#) ElseIf (Alfa < 2# * (Sqr2 - 1#)) Then ' E' data Alfa; calcola la FCt (se < 0.5) corrispondente: Om = Acos(((Alfa ^ 2) + 2# * Alfa - 2#) / (2# * Alfa - 2#)) FCt = Om / PI2 End If A0 = Alfa A1 = 0# B0 = 1# B1 = -(1# - Alfa) End Select ' NCel = 1 NK = 1 ReDim Ac(0 To NK, 1 To NCel) ' Coefficienti del filtro ReDim Bc(0 To NK, 1 To NCel) ' calcolati in Doppia Precisione. ReDim w(0 To NK, 1 To NCel) Ac(0, 1) = A0 Ac(1, 1) = A1 Bc(0, 1) = B0 Bc(1, 1) = B1 AttTot = 1# ' GoTo SintesiIIR_BT_End ' ElseIf (Trans$ = "Resonator") _ Or (Trans$ = "Notch Filter") Then ' Calcolo di tipi di filtro elementari del 2° ordine. ' Vedi Sub zNote: Om = PI2 * FCt ' Pulsazione di risonanza [rad/s]. Cm = Om / Tan(PI * FCt) ' Costante di mappatura. ' ' Calcola i poli: a = 1# + Cm / (q * Om) + (Cm ^ 2 / Om ^ 2) B0 = 1# B1 = 2# * (1# - Cm ^ 2 / Om ^ 2) / a B2 = (1# - Cm / (q * Om) + (Cm ^ 2 / Om ^ 2)) / a ' Select Case Trans$ Case "Resonator" ' Mette due zeri reali a ± 1: A0 = Cm / a A1 = 0# A2 = -Cm / a ' Case "Notch Filter" ' Mette una coppia di zeri complessi coniugati: Rn = Sqr(B2) ' TETAP = Acos(-B1 / (2# * Rn)) A0 = 1# A1 = B1 / Rn ' -2 * Cos(TETAP) A2 = 1# ' Cos(TETAP) ^ 2 + Sin(TETAP) ^ 2 End Select ' NCel = 1 NK = 2 ReDim Ac(0 To NK, 1 To NCel) ' Coefficienti del filtro ReDim Bc(0 To NK, 1 To NCel) ' calcolati in Doppia Precisione. ReDim w(0 To NK, 1 To NCel) Ac(0, 1) = A0 Ac(1, 1) = A1 Ac(2, 1) = A2 Bc(0, 1) = B0 Bc(1, 1) = B1 Bc(2, 1) = B2 AttTot = 1# ' GoTo SintesiIIR_BT_End End If ' '-- Filtri Butterworth e Chebyshev ----------------------------------------------------- FCTA = Tan(PI * FCt) FATA = Tan(PI * FAt) FF = FATA / FCTA If (Proty$ = "Chebyshev") Then Eps = Sqr(1# / (1# - Rip / 100#) ^ 2# - 1#) ATL = 10# ^ (Att / 10#) ' ' COMPUTATION OF THE ORDER OF THE FILTER. ' --------------------------------------- ' N = 0 Do N = N + 1 Select Case Proty$ Case "Butterworth" res# = 1# / (1# + FF ^ (2 * N)) ' Case "Chebyshev" CEB(1) = 1# CEB(2) = FF If (1 < N) Then f = FF + FF For I = 2 To N CEB(I + 1) = f * CEB(I) - CEB#(I - 1) Next I End If res = 1# / (1# + (Eps * CEB(N + 1)) ^ 2) End Select Loop While res >= ATL ' If ((N Mod 2) = 1) Then N = N + 1 Select Case Proty$ Case "Butterworth" res = 1# / (1# + FF ^ (2 * N)) ' Case "Chebyshev" CEB(1) = 1 CEB(2) = FF f = FF + FF For I = 2 To N CEB(I + 1) = f * CEB(I) - CEB(I - 1) Next I res = 1# / (1# + (Eps * CEB(N + 1)) ^ 2) End Select End If ' If (2 * NCMax < N) Then M$ = "**** ERRORE ****" & vbNewLine M$ = M$ & vbNewLine & "Per ottenere il filtro richiesto" & vbNewLine M$ = M$ & "si devono usare " & Str$(N / 2) & " Sezioni." & vbNewLine M$ = M$ & "Il numero Massimo e' " & Str$(NCMax) & "." & vbNewLine M$ = M$ & vbNewLine & "Provare a variare i dati del filtro." MsgBox M$, vbCritical, " modDSP: Sintesi IIR" SintesiIIR_BT = False DoEvents Exit Function End If ' AttTot = 10# * Log(res) / Log10 NCel = N / 2 NK = Switch(Trans$ = "No Transformation", 2, _ Trans$ = "Low Pass -> Low Pass", 2, _ Trans$ = "Low Pass -> High Pass", 2, _ Trans$ = "Low Pass -> Band Pass", 4, _ Trans$ = "Low Pass -> Band Stop", 4) ReDim Ac(0 To NK, 1 To NCel) ' Coefficienti del filtro. ReDim Bc(0 To NK, 1 To NCel) ' Calcolati in Doppia Precisione. ReDim POLRE(1 To NCel) ' Vettori in uso durante la Sintesi. ReDim POLIM(1 To NCel) ' " " " " ReDim w(0 To NK, 1 To NCel) ' ' POLE POSITION EVALUATION. ' ------------------------- ' AT = FCTA Bt = FCTA If (Proty$ = "Chebyshev") Then AB = ((1# + Sqr(1# + Eps ^ 2)) / Eps) ^ (1# / N) AT = (AB - 1# / AB) * AT / 2# Bt = (AB + 1# / AB) * Bt / 2# End If For NC = 1 To NCel TETAP = PI# * (2# * NC - 1#) / (2# * N) POLRE(NC) = -AT * Cos(TETAP) POLIM(NC) = Bt * Sin(TETAP) Next NC For NC = 1 To NCel P1 = POLIM(NC) ^ 2 P2 = (1# - POLRE(NC)) ^ 2 P3 = POLRE(NC) ^ 2 POLRE(NC) = (1# - P3 - P1) / (P2 + P1) POLIM(NC) = (2# * POLIM(NC)) / (P2 + P1) Next NC ' ' EVALUATION OF THE SECOND ORDER SECTION COEFFICIENTS. ' ---------------------------------------------------- ' For NC = 1 To NCel Bc(0, NC) = 1# Bc(1, NC) = -2# * POLRE(NC) Bc(2, NC) = POLRE(NC) * POLRE(NC) + POLIM(NC) * POLIM(NC) Next NC ' Select Case Proty$ Case "Butterworth" COSTA = 1# ' Case "Chebyshev" COSTA = (1# - Rip / 100#) ^ (1# / NCel) End Select ' For NC = 1 To NCel TOT = 4# / (Bc(0, NC) + Bc(1, NC) + Bc(2, NC)) CC = TOT / COSTA Ac(0, NC) = 1# / CC Ac(1, NC) = 2# / CC Ac(2, NC) = 1# / CC Next NC ' ' FREQUENCY TRANSFORMATION. ' ------------------------- ' Select Case Trans$ Case "No Transformation" GoTo SintesiIIR_BT_End ' Case "Low Pass -> Low Pass" AL = Sin(PI * (FCt - F1)) / Sin(PI * (FCt + F1)) A0 = -AL A1 = 1# A2 = 0# B0 = A1 B1 = A0 B2 = 0# ' Case "Low Pass -> High Pass" AL = -Cos(PI * (FCt + F1)) / Cos(PI * (FCt - F1)) A0 = -AL A1 = -1# A2 = 0# B0 = -A1 B1 = -A0 B2 = 0# ' Case "Low Pass -> Band Pass" AL = Cos(PI * (F2 + F1)) / Cos(PI * (F2 - F1)) CAPPA = Tan(PI * FCt) / Tan(PI * (F2 - F1)) A0 = -(CAPPA - 1#) / (CAPPA + 1#) A1 = 2# * AL * CAPPA# / (CAPPA + 1#) A2 = -1# B0 = -A2 B1 = -A1 B2 = -A0 ' Case "Low Pass -> Band Stop" AL = Cos(PI * (F2 + F1)) / Cos(PI * (F2 - F1)) CAPPA = Tan(PI * FCt) * Tan(PI * (F2 - F1)) A0 = (1# - CAPPA) / (1# + CAPPA) A1 = -2# * AL / (1# + CAPPA) A2 = 1# B0 = A2 B1 = A1 B2 = A0 End Select ' AN(0, 0) = B0 * B0 AN(1, 0) = 2# * B0 * B1 AN(2, 0) = B1 * B1 + 2# * B0 * B2 AN(3, 0) = 2# * B1 * B2 AN(4, 0) = B2 * B2 AN(0, 1) = A0 * B0 AN(1, 1) = A0 * B1 + A1 * B0 AN(2, 1) = A0 * B2 + A1 * B1 + A2 * B0 AN(3, 1) = A1 * B2 + A2 * B1 AN(4, 1) = A2 * B2 AN(0, 2) = A0 * A0 AN(1, 2) = 2# * A0 * A1 AN(2, 2) = A1 * A1 + 2# * A0 * A2 AN(3, 2) = 2# * A1 * A2 AN(4, 2) = A2 * A2 ' For NC = 1 To NCel For K = 0 To 2 FINA(K) = Ac(K, NC) FINB(K) = Bc(K, NC) Next K For K = 0 To NK Ac(K, NC) = 0# Bc(K, NC) = 0# For I = 0 To 2 Ac(K, NC) = AN(K, I) * FINA(I) + Ac(K, NC) Bc(K, NC) = AN(K, I) * FINB(I) + Bc(K, NC) Next I Next K Next NC For NC = 1 To NCel Rn = Bc(0, NC) For K = 0 To NK Ac(K, NC) = Ac(K, NC) / Rn Bc(K, NC) = Bc(K, NC) / Rn Next K Next NC ' ' SintesiIIR_BT_End: SintesiIIR_BT_ERR: SintesiIIR_BT = (Err.Number = 0) If (Err.Number <> 0) Then M$ = FiltroErrIIR$ & vbNewLine & vbNewLine M$ = M$ & "Errore " & Str$(Err.Number) & vbNewLine M$ = M$ & Err.Description MsgBox M$, vbCritical, " modDSP: SintesiIIR_BT" End If ' ' ' End Function Private Sub zNote() ' ' Metodi di calcolo dei filtri elementari. ' - Algoritmi di trasformazione da "Digital Signal Processing" ' di W. D. Stanley - pg. 172, 173, 174. ' - La funzione di trasferimento del risuonatore reale parallelo ' e' tratta da "Teoria delle Reti Elettriche", Appunti dai ' corsi del Politecnico di Torino - pg. (1.3)1 e seg. ' - L' idea della sostituzione degli zeri per il "Notch Filter" ' proviene da: http://www-users.cs.york.ac.uk/~fisher/mkfilter/res.html ' '--------------------------------------------------------------------------------------- ' "RC Type Low Pass Filter": ' La funzione di trasferimento e': ' ' U(p) 1 con: ' ---- = ----------- Fc = 1/(2*PI*RC) frequenza di taglio a -3 dB [Hz]. ' I(p) 1 + RC*p C, R capacita' e resistenza ' del circuito. ' ' Le formule di trasformazione (1) da: ' ' A0 + A1*p 1 ' G(p) = --------- --> ------------------- ' B0 + B1*p 1 + (1/(2*PI*Fc))*p ' a: ' -1 ' a0 + a1*z ' H(z) = ----------- ' -1 ' 1 + b1*z ' ' sono: con: ' Cm = 2*PI*FCt/Tan(2*PI*FCt/2) FCt = Fc/Fs --> 0 <= FCt <= 0.5 ' A = B0 + B1*Cm A0 = 1 ' a0 = (A0 + A1*Cm)/A A1 = 0 ' a1 = (A0 - A1*Cm)/A B0 = 1 ' b1 = (B0 - B1*Cm)/A B1 = 1/(2*PI*FCt) ' ' '--------------------------------------------------------------------------------------- ' "RC Type High Pass Filter": ' La funzione di trasferimento e': ' ' U(p) RCp con: ' ---- = ----------- Fc = 1/(2*PI*RC) frequenza di taglio a -3 dB [Hz]. ' I(p) 1 + RC*p C, R capacita' e resistenza ' del circuito. ' ' si usano le formule di trasformazione (1) con: ' FCt = Fc/Fs --> 0 <= FCt <= 0.5 ' A0 = 0 ' A1 = 1/(2*PI*FCt) ' B0 = 1 ' B1 = 1/(2*PI*FCt) ' ' '--------------------------------------------------------------------------------------- ' "Resonator": ' La funzione di trasferimento e': ' ' U(p) P con: ' ---- = --------------------------------- Oc = 1/Sqr(L*C) pulsazione di risonanza [rad/s]. ' I(P) 1 + (1/(Q*Oc))*P + (1/(Oc^2))*P^2 Fc = Oc/(2*PI) frequenza di risonanza [Hz]. ' L, C, Rp induttanza, capacita' e ' resistenza in parallelo ' del circuito risonante. ' Q = Oc*Rp*C coefficiente di risonanza. ' ' Le formule di trasformazione (2) da: ' ' A0 + A1*p + A2*p^2 p ' G(p) = ------------------ --> --------------------------------- ' B0 + B1*p + B2*p^2 1 + (1/(Q*Oc))*P + (1/(Oc^2))*P^2 ' a: ' -1 -2 ' a0 + a1*z + a2*z ' H(z) = ------------------ ' -1 -2 ' 1 + b1*z + b2*z ' ' sono: con: ' Cm = 2*PI*FCt/Tan(2*PI*FCt/2) FCt = Fc/Fs --> 0 <= FCt <= 0.5 ' A = B0 + B1*Cm + B2*Cm^2 A0 = 0 ' a0 = (A0 + A1*Cm + A2*Cm^2)/A A1 = 1 ' a1 = (2*A0 - 2*A2*Cm^2)/A A2 = 0 ' a2 = (A0 - A1*Cm + A2*Cm^2)/A B0 = 1 ' b1 = (2*B0 - 2*B2*Cm^2)/A B1 = 1/(2*PI*FCt*Q) ' b2 = (B0 - B1*Cm + B2*Cm^2)/A B2 = 1/(2*PI*FCt)^2 ' '--------------------------------------------------------------------------------------- ' "Notch Filter": ' L' algoritmo precedente costruisce una funzione di trasferimento H(z) avente ' una coppia di zeri reali a ± 1; questo consente di avere guadagno zero alle ' basse ed alle alte frequenze. ' Il "Notch Filter" viene ricavato dal "Resonator" sostituendo questa coppia di ' zeri reali con due zeri complessi coniugati, giacenti sul cerchio unitario e ' con argomento uguale a quello dei poli (anch' essi complessi coniugati). ' ' Se Zp = r*Exp(± i*Theta) sono i poli di H(z), gli zeri devono essere: ' ' Zz = Exp(± i*Theta) ' ' Viene usata la relazione: ' ' -1 -1 ' (1 - r*Exp(+i*Theta)*z )*(1 - r*Exp(-i*Theta)*z ) = ' ' -1 2 -2 -1 -2 ' = 1 - 2*r*Cos(Theta)*z + r * z = 1 + b1*z + b2*z (i = Sqr(-1)) ' ' da cui: ' ' r = Sqr(b2) ' Theta = Acos(-b1/(2*Sqr(b2)) ' ' e quindi: ' ' -1 -2 -1 -2 ' a0 + a1*z + a2*z = 1 + (b1/(2*Sqr(b2))*z + z ' '--------------------------------------------------------------------------------------- ' "U(t) = a·I(t) + (1-a)·U(t-1)" ' La funzione di trasferimento e': ' ' a ' H(z) = ---------------- ' -1 ' 1 - (1 - a)*z ' ' Usando, invertite, le formule di trasformazione (1) si ha: ' ' a*Cm - a*p ' G(p) = --------------- con: Cm = Oc/Tan(Oc/2) ' a*Cm + (2 - a)*p Oc = 2*PI*Fc ' ' La pulsazione di taglio Oc si ha per: ' ' |a*Cm - a*i*Oc| 1 ' ----------------------- = ------ con: p = i*Oc ' |a*Cm + (2 - a)*i*Oc| Sqr(2) i = Sqr(-1) ' ' da cui si ricava il valore di a necessario per ' ottenere la frequenza di taglio FCt voluta: ' ' a = Cos(Oc) - 1 + Sqr(Cos(Oc)^2 - 4*Cos(Oc) + 3) sempre con: Oc = 2*PI*FCt ' FCt = Fc/Fs ' ' Se invece viene dato a, si ottiene una FCt corrispondente: ' ' FCt = Acos(((a^2) + 2*a - 2)/(2*a - 2))/(2*PI) solo per: a < 2*(Sqr(2) - 1) ' e: FCt < 0.5 ' ' ' End Sub Sub DFT(ByVal NVAL&, d() As Double, r() As Double, X() As Double, _ s() As Double, f() As Double, ByVal NFRE&) ' ' Definizioni: ' ' DFT = Discrete Fourier Transform. ' ' Entra con: ' NVAL = Numero max. dei dati in D() da usare. ' D() = Vettore dei dati Reali del Segnale(t). ' NFRE = Numero dei valori di Frequenza calcolati. ' ' Esce con: ' R() = Vettore dei valori Reali della Trasformata(f). ' X() = Vettore dei valori Immaginari della Trasformata(f). ' S() = Vettore dello Spettro di Amp. della Trasformata(f). ' F() = Vettore delle Fasi della Trasformata(f) [rad]. ' ' I valori della trasformata in R(), X() NON vengono ' normalizzati sulla lunghezza del Segnale NVAL. ' ' Gli NVAL dati nel vettore D() sono organizzati come D(0 To NVAL - 1). ' Gli NFRE dati nei vettori R(), X(), S() e F() sono organizzati come Vettore(0 To NFRE). ' ' Viene calcolata in R(0 To NFRE), X(0 To NFRE) solo la prima ' meta' dello spettro per i valori di frequenza NFRE = NVAL/2 ' piu' la Componente Continua in R(0). ' ' Algoritmi tratti da: "Digital Time Series Analysis" di R. Otnes ' e L. Enochson - pg. 138. Nessun tentativo e' stato fatto per ' migliorare la velocita' di esecuzione (tipo look-up table e/o ' valutazione ricorsiva delle funzioni trigonometriche). ' Dim JF&, JD&, PI2_V#, Omega#, OmT#, NFRE_D# ' 'NFRE = Int(NVAL / 2) NFRE_D = CDbl(NVAL) / 2# PI2_V = PI2 / CDbl(NVAL) ' For JF = 0 To NFRE Omega = PI2_V * CDbl(JF) r(JF) = 0# X(JF) = 0# ' For JD = 0 To NVAL - 1 OmT = Omega * CDbl(JD) r(JF) = r(JF) + d(JD) * Cos(OmT) X(JF) = X(JF) + d(JD) * Sin(OmT) Next JD Next JF ' ' Calcolo dello Spettro di Ampiezza e delle Fasi: s(0) = Abs(r(0)) / CDbl(NVAL) f(0) = IIf(r(0) < 0#, PI, 0#) For JF = 1 To NFRE - 1 s(JF) = Sqr(r(JF) * r(JF) + X(JF) * X(JF)) / NFRE_D f(JF) = DATAN2(-X(JF), r(JF)) ' Scala da -PI a +PI. 'F(JF) = Atan2(-X(JF), R(JF)) ' Scala da 0 a +2PI. Next JF If (NVAL Mod 2) = 0 Then s(NFRE) = Abs(r(NFRE)) / CDbl(NVAL) f(NFRE) = IIf(r(NFRE) < 0#, PI, 0#) Else s(NFRE) = Sqr(r(NFRE) * r(NFRE) + X(NFRE) * X(NFRE)) / NFRE_D f(NFRE) = DATAN2(-X(NFRE), r(NFRE)) ' Scala da -PI a +PI. 'F(NFre) = Atan2(-X(NFre), R(NFre)) ' Scala da 0 a +2PI. End If ' ' ' End Sub Public Sub FFT(d() As Double, r() As Double, X() As Double, _ s() As Double, f() As Double, ByRef NVAL As Long, ByRef NFRE As Long) ' ' Definizioni: ' ' FFT = Fast Fourier Transform. ' ' Entra con: ' D() = Vettore dei dati Reali del Segnale(t). ' NVAL = Numero max. dei dati in D() da usare. ' ' Esce con: ' R() = Vettore dei valori Reali della Trasformata(f). ' X() = Vettore dei valori Immaginari della Trasformata(f). ' S() = Vettore dello Spettro di Amp. della Trasformata(f). ' F() = Vettore delle Fasi della Trasformata(f) [Radianti]. ' NVAL = Numero dei dati usati dalla Trasformata(f). ' NFRE = Numero dei valori di Frequenza calcolati. ' ' I valori della Trasformata in R(), X() NON vengono ' normalizzati sulla lunghezza del Segnale 2^M. ' ' I vettori R(), X(), S() e F() sono dimensionati in questa routine. ' Gli NN dati nei vettori sono organizzati come Vettore(0 To NN - 1). ' ' Tradotta in Basic da: "Theory and Application of Digital ' Signal Processing" di L. Rabiner e B. Gold. - pg. 367. ' ' Ver. 18/10/2005 modificata per AudioCardDSP. ' Dim M&, I&, NMN1&, J&, K&, L&, LE&, LE1&, IP&, JF& Dim TSwap#, Ur#, Ui#, wr#, wi#, Tr#, Ti# ' M = Int(Log(CDbl(NVAL) + 0.5) / Log2) NVAL = 2 ^ M NFRE = NVAL / 2 NMN1 = NVAL - 1 ' ReDim r(0 To NVAL - 1), X(0 To NVAL - 1) ReDim s(0 To NFRE), f(0 To NFRE) ' ' Solo per segnali di ingresso reali e se non si vogliono ' sfruttare le proprieta' "transform-in-place" di questa ' implementazione: For I = 0 To NVAL - 1 r(I) = d(I) X(I) = 0# ' Per segnali di ingresso complessi ' mettere la componente immaginaria. Next I ' J = 0 For I = 0 To NMN1 - 1 If (I < J) Then ' dSwap R(J), R(I): TSwap = r(J) r(J) = r(I) r(I) = TSwap 'dSwap X(J), X(I) ' Solo per segnali complessi. End If K = NFRE Do While (K - 1 < J) J = J - K K = K / 2 Loop J = J + K Next I ' For L = 1 To M LE = 2 ^ L LE1 = LE / 2 Ur = 1# Ui = 0# wr = Cos(PI / LE1) wi = Sin(PI / LE1) For J = 0 To LE1 - 1 For I = J To NVAL - 1 Step LE IP = I + LE1 Tr = r(IP) * Ur - X(IP) * Ui Ti = r(IP) * Ui + X(IP) * Ur r(IP) = r(I) - Tr X(IP) = X(I) - Ti r(I) = r(I) + Tr X(I) = X(I) + Ti Next I Tr = Ur Ti = Ui Ur = Tr * wr - Ti * wi Ui = Tr * wi + Ti * wr Next J Next L ' ' Calcolo dello Spettro di Ampiezza e delle Fasi: s(0) = Abs(r(0)) / NVAL f(0) = IIf(r(0) < 0#, PI, 0#) For JF = 1 To NFRE - 1 s(JF) = Sqr(r(JF) * r(JF) + X(JF) * X(JF)) / NFRE f(JF) = DATAN2(-X(JF), r(JF)) ' Scala da -PI a +PI. 'F(JF) = Atan2(-X(JF), R(JF)) ' Scala da 0 a +2PI. Next JF s(NFRE) = Abs(r(NFRE)) / NVAL f(NFRE) = IIf(r(NFRE) < 0#, PI, 0#) ' ' ' End Sub Public Function WinProf(ByVal WNome$, ByVal NWind&, Optional ByVal Par As Double) As Double() ' ' Calcola e ritorna un vettore con i coefficienti della Window richiesta. ' WNome$: nome della Window richiesta. ' NWind: N° di coefficienti richiesto. ' Par: parametro richiesto per certi tipi di Window. ' Dim K&, KK&, NOdd&, NWind1& Dim BB#, AA#, Ak#, Arg#, ArgL#, ArgW#, NWind_1# ReDim WF_I#(0 To NWind - 1) ' NOdd = NWind Mod 2 NWind1 = Int(NWind / 2) If (NOdd = 1) Then WF_I(NWind1) = 1# NWind_1 = CDbl(NWind - 1) ' AA = 0.5 * (1# - CDbl(NOdd)) BB = CDbl(NWind1) - AA ' ' Calcolo la prima meta' del vettore: Select Case WNome$ Case "Bartlett" ' Zero valued end-points: For K = 0 To NWind1 - 1 WF_I(K) = (2# / NWind_1) * ((NWind_1 / 2#) - Abs(K - (NWind_1 / 2#))) Next K ' Case "Bartlett-Hann" For K = 0 To NWind1 - 1 WF_I(K) = 0.62 - 0.48 * (Abs((K / NWind_1) - 0.5)) _ - 0.38 * Cos(PI2 * CDbl(K) / NWind_1) Next K ' Case "Blackman" For K = 0 To NWind1 - 1 WF_I(K) = 0.42 _ - 0.5 * Cos(PI2 * CDbl(K) / NWind_1) _ + 0.08 * Cos(2# * PI2 * CDbl(K) / NWind_1) Next K ' Case "Blackman-Harris" For K = 0 To NWind1 - 1 WF_I(K) = 0.35875 _ - 0.48829 * Cos(1# * PI2 * CDbl(K) / NWind_1) _ + 0.14128 * Cos(2# * PI2 * CDbl(K) / NWind_1) _ - 0.01168 * Cos(3# * PI2 * CDbl(K) / NWind_1) Next K ' Case "Blackman-Nuttal" For K = 0 To NWind1 - 1 WF_I(K) = 0.3635819 _ - 0.4891775 * Cos(1# * PI2 * CDbl(K) / NWind_1) _ + 0.1365995 * Cos(2# * PI2 * CDbl(K) / NWind_1) _ - 0.0106411 * Cos(3# * PI2 * CDbl(K) / NWind_1) Next K ' Case "Flat top" For K = 0 To NWind1 - 1 WF_I(K) = 0.21557895 _ - 0.41663158 * Cos(1# * PI2 * CDbl(K) / NWind_1) _ + 0.2772631 * Cos(2# * PI2 * CDbl(K) / NWind_1) _ - 0.083578947 * Cos(3# * PI2 * CDbl(K) / NWind_1) _ + 0.006947368 * Cos(4# * PI2 * CDbl(K) / NWind_1) Next K ' Case "Gauss" For K = 0 To NWind1 - 1 WF_I(K) = Exp(-0.5 * ((K - NWind_1 / 2) / (Par * NWind_1 / 2)) ^ 2) Next K ' Case "Hamming generalizzata" For K = 0 To NWind1 - 1 WF_I(K) = Par - (1# - Par) * Cos(PI2 * CDbl(K) / NWind_1) Next K ' Case "Hamming" ' E' l' "Hamming generalizzata" con Par = 0.54: For K = 0 To NWind1 - 1 WF_I(K) = 0.54 - 0.46 * Cos(PI2 * CDbl(K) / NWind_1) Next K ' Case "Hanning" ' E' l' "Hamming generalizzata" con Par = 0.5: For K = 0 To NWind1 - 1 WF_I(K) = 0.5 - 0.5 * Cos(PI2 * CDbl(K) / NWind_1) Next K ' Case "Kaiser" For K = 0 To NWind1 - 1 Ak = Par * Sqr(1# - ((2# * CDbl(K) / NWind_1) - 1#) ^ 2) WF_I(K) = I0_Kaiser(Ak) / I0_Kaiser(Par) Next K ' Case "Lanczos" ArgL = PI / BB For K = 0 To NWind1 - 1 Ak = CDbl(NWind1 - K) - AA Arg = ArgL * Ak WF_I(K) = Abs(Sin(Arg) / Arg) ^ Par Next K ' Case "Nuttal" For K = 0 To NWind1 - 1 WF_I(K) = 0.355768 _ - 0.487396 * Cos(1# * PI2 * CDbl(K) / NWind_1) _ + 0.144232 * Cos(2# * PI2 * CDbl(K) / NWind_1) _ - 0.012604 * Cos(3# * PI2 * CDbl(K) / NWind_1) Next K ' Case "Rettangolare" For K = 0 To NWind1 - 1 WF_I(K) = 1# Next K ' Case "Triangolare" For K = 0 To NWind1 - 1 'WF_I(K) = CDbl(K) / CDbl(NWind1) ' ' Non-zero end-points: WF_I(K) = (2# / NWind) * ((NWind / 2#) - Abs(K - (NWind_1 / 2#))) Next K ' Case "Weber" ArgW = 1.5 * Par / BB For K = 0 To NWind1 - 1 Ak = CDbl(NWind1 - K) - AA Arg = ArgW * Ak If (Arg <= 0.75) Then WF_I(K) = A0 + A1 * Arg + A2 * Arg ^ 2 + A3 * Arg ^ 3 Else WF_I(K) = B0 + B1 * Arg + B2 * Arg ^ 2 + B3 * Arg ^ 3 End If Next K ' Case "Welch" For K = 0 To NWind1 - 1 WF_I(K) = 1# - ((K - NWind_1 / 2#) / (NWind_1 / 2#)) ^ 2 Next K End Select ' ' ' Completo il vettore ' con la seconda meta': For K = 0 To NWind1 - 1 KK = NWind - 1 - K WF_I(KK) = WF_I(K) Next K ' WinProf = WF_I() ' ' ' End Function Private Function I0_Kaiser(ByVal X As Double, Optional ByRef N As Long) As Double ' ' Calcola e ritorna la funzione di Bessel modificata del primo ' tipo e di ordine zero per 0 <= X <= 20. N e' il numero ' di iterazioni impiegate per la convergenza. ' Da: Theory and Application of Digital Signal Processing. ' L. R. Rabiner, B. Gold - pg. 103 ' Dim Y#, e#, de#, SDE# Const t# = 0.00000001 ' Y = X / 2# e = 1# de = 1# For N = 1 To 25 de = de * Y / CDbl(N) SDE = de ^ 2 e = e + SDE If (0# < (e * t - SDE)) Then I0_Kaiser = e Exit Function End If Next N ' N = 0 ' ' ' End Function Public Function WinTipi() As Window_Type() ' ' Imposta i parametri dei tipi di ' "Window" disponibili: ' Profilo da File: WinTipi(0) -> non usato, in questa applicazione. ' Bartlett: WinTipi(1) ' Bartlett-Hann: WinTipi(2) ' Blackman: WinTipi(3) ' Blackman-Harris: WinTipi(4) ' Blackman-Nuttal: WinTipi(5) ' Flat top: WinTipi(6) ' Gauss: WinTipi(7) ' Hamming generalizzata: WinTipi(8) ' Hamming: WinTipi(9) ' Hanning: WinTipi(10) ' Kaiser: WinTipi(11) ' Lanczos: WinTipi(12) ' Nuttal: WinTipi(13) ' Rettangolare: WinTipi(14) ' Triangolare: WinTipi(15) ' Weber: WinTipi(16) ' Welch: WinTipi(17) ' Dim I&, WTipi() As Window_Type ' I = 0 ReDim WTipi(0 To I) WTipi(I).Nome = "--> Nome del File" WTipi(I).PMin = -1 WTipi(I).PMax = -1 ' I = I + 1 ReDim Preserve WTipi(0 To I) WTipi(I).Nome = "Bartlett" WTipi(I).PMin = -1 WTipi(I).PMax = -1 ' I = I + 1 ReDim Preserve WTipi(0 To I) WTipi(I).Nome = "Bartlett-Hann" WTipi(I).PMin = -1 WTipi(I).PMax = -1 ' I = I + 1 ReDim Preserve WTipi(0 To I) WTipi(I).Nome = "Blackman" WTipi(I).PMin = -1 WTipi(I).PMax = -1 ' I = I + 1 ReDim Preserve WTipi(0 To I) WTipi(I).Nome = "Blackman-Harris" WTipi(I).PMin = -1 WTipi(I).PMax = -1 ' I = I + 1 ReDim Preserve WTipi(0 To I) WTipi(I).Nome = "Blackman-Nuttal" WTipi(I).PMin = -1 WTipi(I).PMax = -1 ' I = I + 1 ReDim Preserve WTipi(0 To I) WTipi(I).Nome = "Flat top" WTipi(I).PMin = -1 WTipi(I).PMax = -1 ' I = I + 1 ReDim Preserve WTipi(0 To I) WTipi(I).Nome = "Gauss" WTipi(I).PMin = 0.05 WTipi(I).PMax = 0.5 WTipi(I).PCor = 0.2 ' I = I + 1 ReDim Preserve WTipi(0 To I) WTipi(I).Nome = "Hamming generalizzata" WTipi(I).PMin = 0.5 WTipi(I).PMax = 1# WTipi(I).PCor = 0.5 ' I = I + 1 ReDim Preserve WTipi(0 To I) WTipi(I).Nome = "Hamming" WTipi(I).PMin = -1 WTipi(I).PMax = -1 ' I = I + 1 ReDim Preserve WTipi(0 To I) WTipi(I).Nome = "Hanning" WTipi(I).PMin = -1 WTipi(I).PMax = -1 ' I = I + 1 ReDim Preserve WTipi(0 To I) WTipi(I).Nome = "Kaiser" WTipi(I).PMin = 0.5 WTipi(I).PMax = 20# WTipi(I).PCor = 10# ' I = I + 1 ReDim Preserve WTipi(0 To I) WTipi(I).Nome = "Lanczos" WTipi(I).PMin = 0.5 WTipi(I).PMax = 3.5 WTipi(I).PCor = 1.4 ' I = I + 1 ReDim Preserve WTipi(0 To I) WTipi(I).Nome = "Nuttal" WTipi(I).PMin = -1 WTipi(I).PMax = -1 ' I = I + 1 ReDim Preserve WTipi(0 To I) WTipi(I).Nome = "Rettangolare" WTipi(I).PMin = -1 WTipi(I).PMax = -1 ' I = I + 1 ReDim Preserve WTipi(0 To I) WTipi(I).Nome = "Triangolare" WTipi(I).PMin = -1 WTipi(I).PMax = -1 ' I = I + 1 ReDim Preserve WTipi(0 To I) WTipi(I).Nome = "Weber" WTipi(I).PMin = 0.5 WTipi(I).PMax = 1# WTipi(I).PCor = 1# ' I = I + 1 ReDim Preserve WTipi(0 To I) WTipi(I).Nome = "Welch" WTipi(I).PMin = 0.5 WTipi(I).PMax = -1 WTipi(I).PCor = -1 ' WinTipi = WTipi() ' ' ' End Function Public Sub SFFTBF(ByVal N As Long, d() As Double, _ ByRef Re() As Double, ByRef Im() As Double, ByVal NFRE As Long) ' ' SFFTBF(D, N, MM, S1, C1, S3, C3, ITAB) ' ' A real-valued, in place, split-radix FFT program ' Real input and output in data array D ' Length is N = 2 ** MM ' Decimation-in-time, cos/sin in second loop with table look-up ' Output in order: ' [ Re(0), Re(1), ..., Re(N/2), Im(N/2-1), ..., Im(1) ] ' ' S1 - array of sin() table (length >= N/8-1) ' C1 - array of cos() table (length >= N/8-1) ' S3 - array of sin() table (length >= N/8-1) ' C3 - array of cos() table (length >= N/8-1) ' ITAB - integer bit reversal table (length >= sqrt(2n) ) ' ' The initialization routine SFFTBI must be called prior to calling ' this routine. SFFTBI need not be called again unless N changes. ' ' Original code (RVFFT) written by H.V. Sorensen, ' Rice University, Oct. 1985 ' ' Modifications made by Steve Kifowit, 26 June 1997 ' -- table look-up of sines and cosines ' -- incorporation of bit reversal table ' -- quick return ' ' Tradotta dal FORTRAN e modificata da F. Languasco 25/08/2005. ' ' Entra con: ' N = Numero max. di dati in D() da usare. ' D() = Vettore dei valori Reali del Segnale(t). ' NFRE = Numero dei valori di Frequenza da calcolare. ' ' Esce con: ' Re() = Vettore dei valori Reali della Trasformata(f). ' Im() = Vettore dei valori Immaginari della Trasformata(f). ' ' I valori della trasformata in Re(), Im() NON vengono ' normalizzati sulla lunghezza del Segnale 2^MM. ' ' Gli N dati nel vettore D() sono organizzati come D(0 To N - 1); ' gli NFre_C + 1 dati nei vettori Re() e Im() sono organizzati come ' Vettore(0 To NFre). ' ' Ver: 27/10/2006. ' Dim J&, I&, K&, Ic&, ID&, I0&, I1&, I2&, I3&, I4&, I5&, I6&, I7&, I8& Dim N1&, N2&, N4&, N8&, NN&, It& Dim XT#, R1#, t1#, T2#, T3#, T4#, T5#, T6# Dim CC1#, SS1#, CC3#, SS3#, NFreq# ' If (N = 1) Then Exit Sub ' ' Sposto i dati di ingresso nel vettore D1() per non distruggerli ' e per avere un vettore a base 1 come richiesto da questa routine: CopyMemory D1(1), d(0), 8 * N ' N1 = ITAB(1) For K = 2 To N1 I0 = N1 * ITAB(K) + 1 I = K J = I0 For It = 2 To ITAB(K) + 1 XT = D1(I) D1(I) = D1(J) D1(J) = XT I = I + N1 J = I0 + ITAB(It) Next It Next K ' Ic = 1 ID = 4 70 For I0 = Ic To N Step ID I1 = I0 + 1 R1 = D1(I0) D1(I0) = R1 + D1(I1) D1(I1) = R1 - D1(I1) Next I0 Ic = 2 * ID - 1 ID = 4 * ID If (Ic < N) Then GoTo 70 ' N2 = 2 NN = N / 2 For K = 2 To MM N2 = N2 * 2 N4 = N2 / 4 N8 = N2 / 8 NN = NN / 2 Ic = 0 ID = N2 * 2 40 For I = Ic To N - 1 Step ID I1 = I + 1 I2 = I1 + N4 I3 = I2 + N4 I4 = I3 + N4 t1 = D1(I4) + D1(I3) D1(I4) = D1(I4) - D1(I3) D1(I3) = D1(I1) - t1 D1(I1) = D1(I1) + t1 If (N4 = 1) Then GoTo 38 I1 = I1 + N8 I2 = I2 + N8 I3 = I3 + N8 I4 = I4 + N8 t1 = (D1(I3) + D1(I4)) / Sqr2 T2 = (D1(I3) - D1(I4)) / Sqr2 D1(I4) = D1(I2) - t1 D1(I3) = -D1(I2) - t1 D1(I2) = D1(I1) - T2 D1(I1) = D1(I1) + T2 38 'CONTINUE Next I Ic = 2 * ID - N2 ID = 4 * ID If (Ic < N) Then GoTo 40 For J = 2 To N8 It = (J - 1) * NN CC1 = C1(It) SS1 = S1(It) CC3 = C3(It) SS3 = S3(It) Ic = 0 ID = 2 * N2 36 For I = Ic To N - 1 Step ID I1 = I + J I2 = I1 + N4 I3 = I2 + N4 I4 = I3 + N4 I5 = I + N4 - J + 2 I6 = I5 + N4 I7 = I6 + N4 I8 = I7 + N4 t1 = D1(I3) * CC1 + D1(I7) * SS1 T2 = D1(I7) * CC1 - D1(I3) * SS1 T3 = D1(I4) * CC3 + D1(I8) * SS3 T4 = D1(I8) * CC3 - D1(I4) * SS3 T5 = t1 + T3 T6 = T2 + T4 T3 = t1 - T3 T4 = T2 - T4 T2 = D1(I6) + T6 D1(I3) = T6 - D1(I6) D1(I8) = T2 T2 = D1(I2) - T3 D1(I7) = -D1(I2) - T3 D1(I4) = T2 t1 = D1(I1) + T5 D1(I6) = D1(I1) - T5 D1(I1) = t1 t1 = D1(I5) + T4 D1(I5) = D1(I5) - T4 D1(I2) = t1 Next I Ic = 2 * ID - N2 ID = 4 * ID If (Ic < N) Then GoTo 36 Next J Next K ' ' Sistemo i valori calcolati in D1() ' nei vettori Re(), Im(): 'For I = 0 To NFre ' Re(I) = D1(I + 1) 'Next I CopyMemory Re(0), D1(1), 8 * (NFRE + 1) ' For I = 1 To NFRE - 1 Im(I) = D1(N - I + 1) Next I ' ' ... End of subroutine SFFTBF ... ' End Sub Public Sub SFFTBF_Corr(d() As Double, ByRef Re_C() As Double, ByRef Im_C() As Double) ' ' SFFTBF(D, N_C, MM_C, S1_C, C1_C, S3_C, C3_C, ITAB_C) ' ' A real-valued, in place, split-radix FFT program ' Real input and output in data array D ' Length is N_C = 2 ** MM_C ' Decimation-in-time, cos/sin in second loop with table look-up ' Output in order: ' [ Re_C(0), Re_C(1), ..., Re_C(N_C/2), Im_C(N_C/2-1), ..., Im_C(1) ] ' ' S1_C - array of sin() table (length >= N_C/8-1) ' C1_C - array of cos() table (length >= N_C/8-1) ' S3_C - array of sin() table (length >= N_C/8-1) ' C3_C - array of cos() table (length >= N_C/8-1) ' ITAB_C - integer bit reversal table (length >= sqrt(2n) ) ' ' The initialization routine SFFTBI_Corr must be called prior to calling ' this routine. SFFTBI_Corr need not be called again unless N1_C or ' N2_C change. ' ' Original code (RVFFT) written by H.V. Sorensen, ' Rice University, Oct. 1985 ' ' Modifications made by Steve Kifowit, 26 June 1997 ' -- table look-up of sines and cosines ' -- incorporation of bit reversal table ' -- quick return ' ' Tradotta dal FORTRAN e modificata da F. Languasco 25/08/2005. ' ' Entra con: ' N_C = Numero max. di dati in D() da usare. ' D() = Vettore dei valori Reali del Segnale(t). ' NFre_C = Numero dei valori di Frequenza da calcolare. ' ' Esce con: ' Re_C() = Vettore dei valori Reali della Trasformata(f). ' Im_C() = Vettore dei valori Immaginari della Trasformata(f). ' ' I valori della trasformata in Re_C(), Im_C() NON vengono ' normalizzati sulla lunghezza del Segnale 2^MM_C. ' ' Gli N_C dati nel vettore D() sono organizzati come D(0 To N_C - 1); ' gli NFre_C + 1 dati nei vettori Re_C() e Im_C() sono organizzati come ' Vettore(0 To NFre_C). ' ' Ver: 22/05/2007 per MAutoCorr_FT e MMutuaCorr_FT. ' Dim J&, I&, K&, Ic&, ID&, I0&, I1&, I2&, I3&, I4&, I5&, I6&, I7&, I8& Dim N1&, N2&, N4&, N8&, NN&, It& Dim XT#, R1#, t1#, T2#, T3#, T4#, T5#, T6# Dim CC1#, SS1#, CC3#, SS3#, NFreq# ' If (N_C = 1) Then Exit Sub ' ' Sposto i dati di ingresso nel vettore D1_C() per non distruggerli ' e per avere un vettore a base 1 come richiesto da questa routine: CopyMemory D1_C(1), d(0), 8 * N_C ' N1 = ITAB_C(1) For K = 2 To N1 I0 = N1 * ITAB_C(K) + 1 I = K J = I0 For It = 2 To ITAB_C(K) + 1 XT = D1_C(I) D1_C(I) = D1_C(J) D1_C(J) = XT I = I + N1 J = I0 + ITAB_C(It) Next It Next K ' Ic = 1 ID = 4 70 For I0 = Ic To N_C Step ID I1 = I0 + 1 R1 = D1_C(I0) D1_C(I0) = R1 + D1_C(I1) D1_C(I1) = R1 - D1_C(I1) Next I0 ' Ic = 2 * ID - 1 ID = 4 * ID If (Ic < N_C) Then GoTo 70 ' N2 = 2 NN = N_C / 2 For K = 2 To MM_C N2 = N2 * 2 N4 = N2 / 4 N8 = N2 / 8 NN = NN / 2 Ic = 0 ID = N2 * 2 40 For I = Ic To N_C - 1 Step ID I1 = I + 1 I2 = I1 + N4 I3 = I2 + N4 I4 = I3 + N4 t1 = D1_C(I4) + D1_C(I3) D1_C(I4) = D1_C(I4) - D1_C(I3) D1_C(I3) = D1_C(I1) - t1 D1_C(I1) = D1_C(I1) + t1 If (N4 = 1) Then GoTo 38 I1 = I1 + N8 I2 = I2 + N8 I3 = I3 + N8 I4 = I4 + N8 t1 = (D1_C(I3) + D1_C(I4)) / Sqr2 T2 = (D1_C(I3) - D1_C(I4)) / Sqr2 D1_C(I4) = D1_C(I2) - t1 D1_C(I3) = -D1_C(I2) - t1 D1_C(I2) = D1_C(I1) - T2 D1_C(I1) = D1_C(I1) + T2 38 'CONTINUE Next I ' Ic = 2 * ID - N2 ID = 4 * ID If (Ic < N_C) Then GoTo 40 ' For J = 2 To N8 It = (J - 1) * NN CC1 = C1_C(It) SS1 = S1_C(It) CC3 = C3_C(It) SS3 = S3_C(It) Ic = 0 ID = 2 * N2 36 For I = Ic To N_C - 1 Step ID I1 = I + J I2 = I1 + N4 I3 = I2 + N4 I4 = I3 + N4 I5 = I + N4 - J + 2 I6 = I5 + N4 I7 = I6 + N4 I8 = I7 + N4 t1 = D1_C(I3) * CC1 + D1_C(I7) * SS1 T2 = D1_C(I7) * CC1 - D1_C(I3) * SS1 T3 = D1_C(I4) * CC3 + D1_C(I8) * SS3 T4 = D1_C(I8) * CC3 - D1_C(I4) * SS3 T5 = t1 + T3 T6 = T2 + T4 T3 = t1 - T3 T4 = T2 - T4 T2 = D1_C(I6) + T6 D1_C(I3) = T6 - D1_C(I6) D1_C(I8) = T2 T2 = D1_C(I2) - T3 D1_C(I7) = -D1_C(I2) - T3 D1_C(I4) = T2 t1 = D1_C(I1) + T5 D1_C(I6) = D1_C(I1) - T5 D1_C(I1) = t1 t1 = D1_C(I5) + T4 D1_C(I5) = D1_C(I5) - T4 D1_C(I2) = t1 Next I ' Ic = 2 * ID - N2 ID = 4 * ID If (Ic < N_C) Then GoTo 36 Next J Next K ' ' Sistemo i valori calcolati in D1_C() ' nei vettori Re_C(), Im_C(): ' For I = 0 To NFre_C - 1 ' Re_C(I) = D1_C(I + 1) ' Next I CopyMemory Re_C(0), D1_C(1), 8 * (NFre_C + 1) ' For I = 1 To NFre_C - 1 Im_C(I) = D1_C(N_C - I + 1) Next I ' ' ... End of subroutine SFFTBF ... ' End Sub Public Sub SFFTBI(ByVal N As Long) ' ' SFFTBI( N, MM, S1, C1, S3, C3, ITAB ) ' ' Table initialization routine for SFFTBF and SFFTBB ' ' Usage: CALL SFFTBI( N, MM, S1, C1, S3, C3, ITAB ) ' Parameters: ' N - integer length of transform (must be a power of two) ' MM - integer such that N = 2**MM ' S1 - array of sin() table (length >= n/8-1) ' C1 - array of cos() table (length >= n/8-1) ' S3 - array of sin() table (length >= n/8-1) ' C3 - array of cos() table (length >= n/8-1) ' ITAB - integer bit reversal table (length >= sqrt(2n)) ' ' Uses standard FORTRAN programs - sin, cos ' ' Steve Kifowit, 26 June 1997 ' ' Tradotta dal FORTRAN e modificata da F. Languasco 25/08/2005. ' ' I valori calcolati nei vettori S1(), C1(), S3(), C3() e ITAB() ed MM ' non vengono ritornati da questa routine ma memorizzati localmente ' per essere, successivamente, usati da SFFTBF. ' ' ' Dim I&, K&, IMax&, M2&, MS& Dim ANG#, t#, u# ' ReDim S1(1 To N / 8 - 1), C1(1 To N / 8 - 1) ReDim S3(1 To N / 8 - 1), C3(1 To N / 8 - 1) ReDim ITAB(1 To Sqr(2 * N)) ReDim D1(1 To N) ' MM = CLng(Log(CDbl(N)) / Log2) ' ' ... Compute bit reversal table ... M2 = Int(MM / 2) MS = 2 ^ M2 If (2 * M2 <> MM) Then M2 = M2 + 1 ITAB(1) = 0 ITAB(2) = 1 IMax = 1 For K = 2 To M2 IMax = 2 * IMax For I = 1 To IMax ITAB(I) = 2 * ITAB(I) ITAB(I + IMax) = 1 + ITAB(I) Next I Next K ITAB(1) = MS ' ' ... Quick return ... If (N <= 8) Then Exit Sub ' ' ... Compute trig tables ... ANG = PI2 / CDbl(N) K = N / 8 - 1 For I = 1 To K t = ANG * CDbl(I) C1(I) = Cos(t) S1(I) = Sin(t) u = 3# * t C3(I) = Cos(u) S3(I) = Sin(u) Next I ' ' ... End of subroutine SFFTBI ... ' End Sub Public Sub SFFTBI_Corr(ByVal N1 As Long, Optional ByVal N2 As Long = -1) ' ' SFFTBI( N, MM_C, S1_C, C1_C, S3_C, C3_C, ITAB_C ) ' ' Table initialization routine for SFFTBF_Corr and SFFTBB_Corr ' ' Usage: CALL SFFTBI( N, MM_C, S1_C, C1_C, S3_C, C3_C, ITAB_C ) ' Parameters: ' N1 - integer length of 1° data vector. ' N2 - integer length of 2° data vector. ' N_C - integer length of transform (is a power of two). ' MM_C - integer such that N_C = 2**MM_C. ' S1_C - array of sin() table (length >= n/8-1). ' C1_C - array of cos() table (length >= n/8-1). ' S3_C - array of sin() table (length >= n/8-1). ' C3_C - array of cos() table (length >= n/8-1). ' ITAB_C - integer bit reversal table (length >= sqrt(2n)). ' ' Uses standard FORTRAN programs - sin, cos ' ' Steve Kifowit, 26 June 1997 ' ' Tradotta dal FORTRAN e modificata da F. Languasco 25/08/2005. ' ' I valori calcolati nei vettori S1_C(), C1_C(), S3_C(), C3_C(), ITAB_C() ' e nelle variabili N_C, MM_C non vengono ritornati da questa routine ma ' memorizzati localmente per essere, successivamente, usati da SFFTBF_Corr ' e SFFTBB_Corr. ' ' Versione per AudioCardDSP (22/05/2007) da usare con MAutoCorr_FT ' e MMutuaCorr_FT. ' Dim I&, K&, IMax&, M2&, MS& Dim OmSh#, ANG#, t#, u# ' N1_C = N1 If (N2 = -1) Then N2_C = N1 Else N2_C = N2 End If ' MM_C = CLng(Ceil(Log(CDbl(N1_C + N2_C - 1)) / Log2)) N_C = 2 ^ MM_C NFre_C = N_C / 2 ' ReDim S1_C(1 To N_C / 8 - 1), C1_C(1 To N_C / 8 - 1) ReDim S3_C(1 To N_C / 8 - 1), C3_C(1 To N_C / 8 - 1) ReDim ITAB_C(1 To Sqr(2 * N_C)) ReDim D1_C(1 To N_C) ReDim Re1_C(0 To NFre_C), Im1_C(0 To NFre_C) ReDim Re2_C(0 To NFre_C), Im2_C(0 To NFre_C) ReDim s(0 To NFre_C), f(0 To NFre_C) ReDim WnRe(0 To NFre_C), WnIm(0 To NFre_C) ' ' Calcolo le tavole dei seni/coseni per ' la trasformata del segnale ritardato ' di N1_C campioni (solo per AutoCorr_FT): OmSh = PI2 * N1_C / N_C For I = 0 To NFre_C WnRe(I) = Cos(OmSh * I) WnIm(I) = -Sin(OmSh * I) Next I ' ' ... Compute bit reversal table ... M2 = Int(MM_C / 2) MS = 2 ^ M2 If (2 * M2 <> MM_C) Then M2 = M2 + 1 ITAB_C(1) = 0 ITAB_C(2) = 1 IMax = 1 For K = 2 To M2 IMax = 2 * IMax For I = 1 To IMax ITAB_C(I) = 2 * ITAB_C(I) ITAB_C(I + IMax) = 1 + ITAB_C(I) Next I Next K ITAB_C(1) = MS ' ' ... Quick return ... If (N_C <= 8) Then Exit Sub ' ' ... Compute trig tables ... ANG = PI2 / CDbl(N_C) K = N_C / 8 - 1 For I = 1 To K t = ANG * CDbl(I) C1_C(I) = Cos(t) S1_C(I) = Sin(t) u = 3# * t C3_C(I) = Cos(u) S3_C(I) = Sin(u) Next I ' ' ... End of subroutine SFFTBI ... ' End Sub Public Sub LeggiCoefficienti_IIR(ByRef NK_U As Long, ByRef NCel_U As Long, _ ByRef W_U() As Single, ByRef Ac_U() As Single, ByRef Bc_U() As Single) ' ' Con questa routine e' possibile ritornare, per uso locale, i ' coefficienti di un filtro IIR calcolati precedentemente dalla ' Function SintesiIIR_BT(): ' ' Routine usata in AudioCardDSP (27/04/2009) ' NK_U = NK ' Ordine NCel_U = NCel ' e numero di sezioni del filtro. W_U() = w() ' Registri delle sezioni del filtro. Ac_U() = Ac() ' Coefficienti del filtro. Bc_U() = Bc() ' " " " ' ' ' End Sub Public Sub LeggiCoefficienti_SPM_R(ByRef Sn_1_U As Complex, ByRef SmRE_U() As Double, _ ByRef z1_1_U As Complex, ByRef z1_N_U As Complex, ByRef NormSpm_U As Double) ' ' Con questa routine e' possibile ritornare, per uso locale, i ' coefficienti di un filtro SPM_R calcolati precedentemente dalla ' Sub SPM_MF_R_Init(): ' ' Routine usata in AudioCardDSP (27/04/2009) ' Sn_1_U = Sn_1 ' Sn * Z^-1 SmRE_U() = SmRE() ' Registro a scorrimento dei campioni del segnale. z1_1_U = z1_1 ' z1(0)^-1 z1_N_U = z1_N ' z1(0)^-N NormSpm_U = NormSpm ' Fattore di normalizzazione sul N° di campioni. ' ' ' End Sub
visual basic sux. learn c++ instead
ReplyDeleteAmazing article; I liked how you talked about the development of top developers and how they are very popular. Finding work is not that tough for them; I also came across a freelancing platform Eiliana.com that provides you with good projects; you can also look at it.
ReplyDelete