Visual Basic 6.0 - Superior Source Code: 2021

Wednesday, October 20, 2021

Ant colony simulation (by David Rutten)

Simulate the behaviour of an Ant colony. This is an experiment of a system with 'intelligent' agents who affect and are affected by their environment and counterparts.

Download from ME
Download from PSC

Sources:
https://planet-source-code.com/vb/scripts/ShowCode.asp?txtCodeId=55195&lngWId=1

Wednesday, September 29, 2021

Voronoi diagrams (by David Rutten)

Voronoi diagrams are named after Russian mathematician Georgy Feodosievych Voronoy who defined and studied the general n-dimensional case in 1908.

Voronoi diagrams are quite useful tools in computational geometry and have a wide range of uses such as, calculating the area per tree in the forest, or figuring out where the poisoned wells were in a city (based on victims' addresses), and so on. In general it is useful for finding "who is closest to whom." A collection of problems where Voronoi diagrams are used is shown below:

Collision detection
Pattern recognition
Geographical optimization
Geometric clustering
Closest pairs algorithms
k-nearest-neighbor queries

Download from ME
Download from PSC

Sources:

https://planet-source-code.com/vb/scripts/ShowCode.asp?txtCodeId=57062&lngWId=1

https://www.codeproject.com/Articles/411242/Create-a-Voronoi-diagram-1-of-2

Saturday, July 31, 2021

CryptoAPI encryption/Decryption, hashing, and random number generation (by Kenneth Ives)

Demonstration of CryptoAPI encryption/Decryption, hashing, and random number generation. Encryption includes RC2, RC4, DES, DES3, AES-128, AES-192, AES-256. Hashing includes MD4, MD5, SHA-1, SHA-256, SHA-384, SHA-512. Windows XP SP2 or earlier no longer supported. Updated 20-Jul-2017, support modules and documentation.

Download from ME
Download from PSC

Count Lines of Code v2.0.353 22 Mar 2017 06:29 AM

CryptoAPI_Group.vbg
Kenneth Ives kenaso@tx.rr.com
------------------------------------------------------------------------------

VBG Name: C:\Kens Software\CryptoAPI\CryptoAPI_Group.vbg

VBP Name: C:\Kens Software\CryptoAPI\CryptoAPIDemo.vbp

Module Name: C:\Kens Software\CryptoAPI\clsKeyEdit.cls

10 Sub routines
80 Miscellaneous lines of code
15 Auto generated lines (** Not included in totals **)
60 Blank lines (** Not included in totals **)
287 Comment lines (** Not included in totals **)
----------
90 Module lines of code

Module Name: C:\Kens Software\CryptoAPI\clsManifest.cls

3 Constant variables
1 Enum Structures
1 Type Structures
5 API Declare statements
4 Property Let routines
2 Sub routines
4 Functions
182 Miscellaneous lines of code
13 Auto generated lines (** Not included in totals **)
88 Blank lines (** Not included in totals **)
289 Comment lines (** Not included in totals **)
----------
198 Module lines of code

Module Name: C:\Kens Software\CryptoAPI\clsOperSystem.cls

105 Constant variables
1 Type Structures
6 API Declare statements
163 Property Get routines
13 Sub routines
1,170 Miscellaneous lines of code
15 Auto generated lines (** Not included in totals **)
381 Blank lines (** Not included in totals **)
378 Comment lines (** Not included in totals **)
----------
1,295 Module lines of code

Module Name: C:\Kens Software\CryptoAPI\clsPrivileges.cls

11 Constant variables
3 Type Structures
6 API Declare statements
4 Property Let routines
1 Property Get routines
3 Sub routines
3 Functions
173 Miscellaneous lines of code
13 Auto generated lines (** Not included in totals **)
117 Blank lines (** Not included in totals **)
213 Comment lines (** Not included in totals **)
----------
199 Module lines of code

Module Name: C:\Kens Software\CryptoAPI\frmAbout.frm

7 Constant variables
3 API Declare statements
9 Sub routines
95 Miscellaneous lines of code
332 Auto generated lines (** Not included in totals **)
23 Blank lines (** Not included in totals **)
84 Comment lines (** Not included in totals **)
----------
114 Module lines of code

Module Name: C:\Kens Software\CryptoAPI\frmMain.frm

6 Constant variables
2 API Declare statements
37 Sub routines
1,340 Miscellaneous lines of code
779 Auto generated lines (** Not included in totals **)
93 Blank lines (** Not included in totals **)
258 Comment lines (** Not included in totals **)
----------
1,385 Module lines of code

Module Name: C:\Kens Software\CryptoAPI\frmSplash.frm

2 Sub routines
10 Miscellaneous lines of code
52 Auto generated lines (** Not included in totals **)
6 Blank lines (** Not included in totals **)
17 Comment lines (** Not included in totals **)
----------
12 Module lines of code

Module Name: C:\Kens Software\CryptoAPI\modCentering.bas

12 Constant variables
1 Type Structures
7 API Declare statements
2 Sub routines
3 Functions
83 Miscellaneous lines of code
1 Auto generated lines (** Not included in totals **)
73 Blank lines (** Not included in totals **)
207 Comment lines (** Not included in totals **)
----------
108 Module lines of code

Module Name: C:\Kens Software\CryptoAPI\modCommon.bas

5 Constant variables
5 Functions
79 Miscellaneous lines of code
1 Auto generated lines (** Not included in totals **)
50 Blank lines (** Not included in totals **)
99 Comment lines (** Not included in totals **)
----------
89 Module lines of code

Module Name: C:\Kens Software\CryptoAPI\modDialogBox.bas

61 Constant variables
7 Type Structures
24 API Declare statements
6 Sub routines
11 Functions
652 Miscellaneous lines of code
1 Auto generated lines (** Not included in totals **)
105 Blank lines (** Not included in totals **)
770 Comment lines (** Not included in totals **)
----------
761 Module lines of code

Module Name: C:\Kens Software\CryptoAPI\modMain.bas

34 Constant variables
1 Type Structures
23 API Declare statements
8 Sub routines
11 Functions
340 Miscellaneous lines of code
1 Auto generated lines (** Not included in totals **)
244 Blank lines (** Not included in totals **)
654 Comment lines (** Not included in totals **)
----------
417 Module lines of code

Module Name: C:\Kens Software\CryptoAPI\modMessages.bas

24 Constant variables
2 Enum Structures
1 Type Structures
9 API Declare statements
4 Sub routines
4 Functions
101 Miscellaneous lines of code
1 Auto generated lines (** Not included in totals **)
79 Blank lines (** Not included in totals **)
405 Comment lines (** Not included in totals **)
----------
145 Module lines of code

Module Name: C:\Kens Software\CryptoAPI\modProcesses.bas

9 Constant variables
1 Type Structures
16 API Declare statements
2 Sub routines
13 Functions
265 Miscellaneous lines of code
1 Auto generated lines (** Not included in totals **)
199 Blank lines (** Not included in totals **)
398 Comment lines (** Not included in totals **)
----------
306 Module lines of code

Module Name: C:\Kens Software\CryptoAPI\modTrimStr.bas

1 Constant variables
1 Type Structures
2 API Declare statements
2 Functions
53 Miscellaneous lines of code
1 Auto generated lines (** Not included in totals **)
33 Blank lines (** Not included in totals **)
129 Comment lines (** Not included in totals **)
----------
59 Module lines of code

----------
5,178 Sub-total for project

VBP Name: C:\Kens Software\CryptoAPI\DLL\kiCryptoAPI.vbp

Module Name: C:\Kens Software\CryptoAPI\DLL\clsAPI_Hash.cls

39 Constant variables
1 Enum Structures
7 API Declare statements
5 Property Let routines
2 Property Get routines
1 Event routines
6 Sub routines
7 Functions
320 Miscellaneous lines of code
15 Auto generated lines (** Not included in totals **)
204 Blank lines (** Not included in totals **)
554 Comment lines (** Not included in totals **)
----------
381 Module lines of code

Module Name: C:\Kens Software\CryptoAPI\DLL\clsBigFiles.cls

19 Constant variables
10 API Declare statements
2 Property Let routines
2 Event routines
4 Sub routines
13 Functions
351 Miscellaneous lines of code
15 Auto generated lines (** Not included in totals **)
230 Blank lines (** Not included in totals **)
704 Comment lines (** Not included in totals **)
----------
399 Module lines of code

Module Name: C:\Kens Software\CryptoAPI\DLL\clsCRC32.cls

8 Constant variables
1 Property Let routines
2 Property Get routines
1 Event routines
3 Sub routines
4 Functions
202 Miscellaneous lines of code
15 Auto generated lines (** Not included in totals **)
56 Blank lines (** Not included in totals **)
189 Comment lines (** Not included in totals **)
----------
218 Module lines of code

Module Name: C:\Kens Software\CryptoAPI\DLL\clsCipher.cls

7 Constant variables
7 Property Let routines
5 Property Get routines
1 Event routines
3 Sub routines
4 Functions
324 Miscellaneous lines of code
15 Auto generated lines (** Not included in totals **)
63 Blank lines (** Not included in totals **)
98 Comment lines (** Not included in totals **)
----------
339 Module lines of code

Module Name: C:\Kens Software\CryptoAPI\DLL\clsCryptoAPI.cls

73 Constant variables
2 Enum Structures
13 API Declare statements
7 Property Let routines
3 Property Get routines
2 Event routines
5 Sub routines
11 Functions
708 Miscellaneous lines of code
15 Auto generated lines (** Not included in totals **)
376 Blank lines (** Not included in totals **)
654 Comment lines (** Not included in totals **)
----------
814 Module lines of code

Module Name: C:\Kens Software\CryptoAPI\DLL\clsHash.cls

1 Constant variables
4 Property Let routines
2 Property Get routines
1 Event routines
3 Sub routines
2 Functions
128 Miscellaneous lines of code
15 Auto generated lines (** Not included in totals **)
38 Blank lines (** Not included in totals **)
117 Comment lines (** Not included in totals **)
----------
135 Module lines of code

Module Name: C:\Kens Software\CryptoAPI\DLL\clsRandom.cls

71 Constant variables
3 Enum Structures
14 API Declare statements
2 Property Let routines
2 Property Get routines
13 Sub routines
19 Functions
918 Miscellaneous lines of code
15 Auto generated lines (** Not included in totals **)
502 Blank lines (** Not included in totals **)
1,391 Comment lines (** Not included in totals **)
----------
1,038 Module lines of code

Module Name: C:\Kens Software\CryptoAPI\DLL\modCommon.bas

18 Constant variables
7 API Declare statements
3 Sub routines
18 Functions
232 Miscellaneous lines of code
1 Auto generated lines (** Not included in totals **)
172 Blank lines (** Not included in totals **)
544 Comment lines (** Not included in totals **)
----------
278 Module lines of code

Module Name: C:\Kens Software\CryptoAPI\DLL\modMessages.bas

24 Constant variables
2 Enum Structures
1 Type Structures
9 API Declare statements
4 Sub routines
4 Functions
101 Miscellaneous lines of code
1 Auto generated lines (** Not included in totals **)
79 Blank lines (** Not included in totals **)
405 Comment lines (** Not included in totals **)
----------
145 Module lines of code

Module Name: C:\Kens Software\CryptoAPI\DLL\modTrimStr.bas

1 Constant variables
1 Type Structures
2 API Declare statements
2 Functions
53 Miscellaneous lines of code
1 Auto generated lines (** Not included in totals **)
33 Blank lines (** Not included in totals **)
129 Comment lines (** Not included in totals **)
----------
59 Module lines of code

----------
3,806 Sub-total for project

==========
8,984 Total number of lines of code

******************************************************************************
NOTE: Visual Basic trailers are not counted. These are the
logical ending statements used by proceedural headings.

End Sub End Function End Property
End If End Type Loop
Next Wend End With
End Select
******************************************************************************

kiCryptoAPI.dll Kenneth Ives (kenaso|at|tx.rr.com)
I am open to ways to improve this application, please email me.
Visual Basic 6.0 with Service Pack 6 runtime files required.
To obtain required files (VBRun60sp6.exe):
http://www.microsoft.com/downloads/details.aspx?FamilyId=7B9BA261-7A9C-43E7-9117-F6730
77FFB3C
VBRun60sp6.exe installs Visual Basic 6.0 SP6 run-time files.
http://support.microsoft.com/kb/290887
This software has been tested on Windows XP SP3 64-bit through Windows 10.
Windows XP 32-bit, 9x, 2000 and NT4 are no longer supported.
*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*
*** WARNING *** WARNING *** WARNING *** WARNING *** WARNING *** WARNING ***
*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*
You acknowledge that this software is subject to the export control
laws and regulations of the United States ("U.S.") and agree to abide
by those laws and regulations. Under U.S. law, this software may not
be downloaded or otherwise exported, reexported, or transferred to
restricted countries, restricted end-users, or for restricted
end-uses. The U.S. currently has embargo restrictions against Cuba,
Iran, Iraq, Libya, North Korea, Sudan, and Syria. The lists of
restricted end-users are maintained on the U.S. Commerce Department's
Denied Persons List, the Commerce Department's Entity List, the
Commerce Department's List of Unverified Persons, and the U.S.
Treasury Department's List of Specially Designated Nationals and
Blocked Persons. In addition, this software may not be downloaded or
otherwise exported, reexported, or transferred to an end-user engaged
in activities related to weapons of mass destruction.
*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*
REFERENCE:
The Cryptography API, or How to Keep a Secret
http://msdn.microsoft.com/en-us/library/ms867086.aspx
CryptoAPI Cryptographic Service Providers
http://msdn.microsoft.com/en-us/library/bb931357(VS.85).aspx
SHA-2 support on MS Windows
Paraphrasing: Regarding SHA-224 support, SHA-224 offers less security
than SHA-256 but takes the same amount of resources. Also SHA-224 is
not generally used by protocols and applications. The NSA's (National
Security Agency) Suite B standards also does not include it. Microsoft
has no plans to add it to future versions of their Cryptographic
Service Providers (CSP).
http://blogs.msdn.com/b/alejacma/archive/2009/01/23/sha-2-support-on-windows-xp.aspx
NIST (National Institute of Standards and Technology)
FIPS (Federal Information Processing Standards Publication)
SP (Special Publications)
http://csrc.nist.gov/publications/PubsFIPS.html
FIPS 180-2 (Federal Information Processing Standards Publication)
dated 1-Aug-2002, with Change Notice 1, dated 25-Feb-2004
http://csrc.nist.gov/publications/fips/fips180-2/FIPS180-2_changenotice.pdf
FIPS 180-3 (Federal Information Processing Standards Publication)

dated Oct-2008 (supercedes FIPS 180-2)
http://csrc.nist.gov/publications/fips/fips180-3/fips180-3_final.pdf
FIPS 180-4 (Federal Information Processing Standards Publication)
dated Mar-2012 (Supercedes FIPS-180-3)
http://csrc.nist.gov/publications/fips/fips180-4/fips-180-4.pdf
Examples of the implementation of the secure hash algorithms
SHA-1, SHA-224, SHA-256, SHA-384, SHA-512, SHA-512/224 and
SHA-512/256, can be found at:
http://csrc.nist.gov/groups/ST/toolkit/examples.html
http://csrc.nist.gov/groups/ST/toolkit/documents/Examples/SHA2_Additional.pdf
Aaron Gifford's additional test vectors
http://www.adg.us/computers/sha.html
Guidelines for Media Sanitization (SP800-88)
http://csrc.nist.gov/publications/nistpubs/800-88/NISTSP800-88_rev1.pdf
WARNING:
MD4 Message-Digest Algorithm has been compromised at the rump
session of Crypto 2004 it was announced that Xiaoyun Wang,
Dengguo Feng, Xuejia Lai and Hongbo Yu found collisions for
MD4, MD5, RIPEMD, and the 128-bit version of HAVAL.
http://eprint.iacr.org/2004/199.pdf
Feb-2005: SHA-1 has been compromised. Recommended that
you do not use for password or document authentication.
http://www.schneier.com/blog/archives/2005/02/sha1_broken.html
http://csrc.nist.gov/groups/ST/toolkit/documents/shs/NISTHashComments-final.pdf
Mar-2005 Demonstrating a technique for finding MD5 collisions quickly.
Eight hours on 1.6 GHz computer.
http://cryptography.hyperlink.cz/md5/MD5_collisions.pdf
Jun-2005 Two researchers from the Institute for Cryptology and
IT-Security have generated PostScript files with identical MD5-sums
but entirely different (but meaningful!) content.
http://www.schneier.com/blog/archives/2005/06/more_md5_collis.html
March 15, 2006: The SHA-2 family of hash functions
(i.e., SHA-224, SHA-256, SHA-384 and SHA-512) may be used
by Federal agencies for all applications using secure hash
algorithms. Federal agencies should stop using SHA-1 for
digital signatures, digital time stamping and other
applications that require collision resistance as soon as
practical, and must use the SHA-2 family of hash functions
for these applications after 2010. After 2010, Federal
agencies may use SHA-1 only for the following applications:
- hash-based message authentication codes (HMACs)
- key derivation functions (KDFs)
- random number generators (RNGs)
Regardless of use, NIST encourages application and protocol
designers to use the SHA-2 family of hash functions for all
new applications and protocols.
http://csrc.nist.gov/groups/ST/hash/policy.html

Export Control: Certain cryptographic devices and technical
data regarding them are subject to Federal export controls.
Exports of cryptographic modules implementing this standard
and technical data regarding them must comply with these
Federal regulations and be licensed by the Bureau of Export
Administration of the U.S. Department of Commerce.

Information about export regulations is available at:
http://www.bis.doc.gov/index.htm
*****************************************************************************
How to use:
For a simple example, execute the SHA_Demo application. The demo converts
the data to a byte array prior to passing it to the DLL to be hashed.
[STRING DATA]
Convert string data to byte array prior to passing to the HashString function.
Ex: abytData() = StrConv("abc", vbFromUnicode)
[FILE DATA]
Just the path and filename are passed in the byte array. Convert the
path\filename data to byte array prior to passing to the HashFile function.
The HashFile routine will open and read the file into an internal byte array.
Ex: abytData() = StrConv("C:\Files\Test Folder\Testfile.txt", vbFromUnicode)
Both will create a hashed output string based on file data input.
-------------------------------------------------------------------------------
Test data provided to test either hash or cipher:
TestPhrase.txt ASCII text phrase (Copy & paste phrase for string test)
TestFile.txt ASCII text file
Binary test files:
kB_32.dat 32,768 binary zeros
OneMil_0.dat One million binary zeros (FIPS 180-3)
OneMil_a.dat One million letter "a" (FIPS 180-2)
API32.txt Text file over 1 MB
*****************************************************************************
Note from Mark Hutchinson's presentation about Microsoft's VB random number
generator. http://www.15seconds.com/issue/051110.htm
References:
Randomize Statement Doesn't Re-initialize Rnd Function
http://support.microsoft.com/default.aspx?scid=kb;en-us;120587
"To re-initialize the random-number generator, use the Rnd function with a
value of -1 to re-initialize the Rnd function, and then use the Randomize
statement with the value you want to use as the seed value for the Rnd
function."

VBA's Pseudo Random Number Generator
http://www.noesis.net.au/prng.php
INFO: How Visual Basic Generates Pseudo-Random Numbers for the RND Function
http://support.microsoft.com/kb/231847/en-us
RND and RANDOMIZE Alternatives for Generating Random Numbers
http://support.microsoft.com/kb/28150/EN-US/

** Enhanced ciphers
********************************
With all ciphers, except ArcFour, the data length will change. After
encrytption, data sizes will not match original sizes. This is due to
internal padding and the storing of information required to decrypt the
data later.
********************************
** PASSWORDS
********************************
Currently there is a minimum and maximum length of the password the user
may enter. This can be changed in the kiCrypt DLL basCommon.bas module.
In the declarations section, locate these two constants and make the
desired change. Be sure to recompile the DLL and the demo application.
PWD_LENGTH_MIN = 8
PWD_LENGTH_MAX = 50
If no hash algorithm is selected then the default hash will be SHA-256.

Sources:
http://www.planetsourcecode.com/vb/scripts/ShowCode.asp?txtCodeId=74645&lngWId=1
http://www.vbforums.com/showthread.php?831741-How-do-I-use-Crypto-API-functionality-in-VB6-without-CAPICOM-ActiveX-control

Wednesday, July 21, 2021

ProcCounters & ProcMonitor - instrument your application

An issue many of us deal with is trying different approaches in an application in order to improve performance. Often enough these are big changes, for example ADO Client vs. Server cursor location or using a file vs. keeping everything in memory.

Download from VBForums
Download from ME

So you don't need a code profiler right off the bat to micro-optimize, instead you need more "global" performance numbers: accumulating counters for the process. There are some API calls to retrieve a number of statistics. Some of the more useful ones measure CPU use, I/O use, and memory use.

ProcCounters is a VB6 class wrapping several of these calls. ProcMonitor is a VB6 UserControl that displays summary information you can watch while running your program. It samples statistics via ProcCounters and shows them in abbreviated format. The test program in the attachment just does a bunch of grinding away while it logs ProcCounters results and has a ProcMonitor (blue here) running as well. These require Windows 2000 or later. Note that ProcMonitor uses SHLWAPI calls to format byte-count values in "base 2" scales, i.e. 1KB = 1024 bytes, etc.

Source:
http://earlier189.rssing.com/browser.php?indx=6373759&item=371

Thursday, July 15, 2021

Digital Signal Processing (VB6)

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