Imports FFTWSharp
Imports System.buffer
Imports System.Runtime.InteropServices
Imports System.math
Imports System.Management
Imports System.Threading
Imports System.Threading.Thread
Public Class SpectralProcessing
Public bytespersample As Integer = 4
Public N As Integer = 256 ' frame size
Public k As Integer = 1 ' number of frames in filter
Public threadcount As Integer ' number of threads
Public corecount As Integer ' number of processor cores
Dim filtertaps()() As Double ' filter taps
Dim currentthread As Integer ' pointer to a thread
Public filterpath As String = ""
Dim integrationtime As Double ' integration time (to set time limit)
' allocate threads
Dim buf2()() As Short ' interpreted as short
Dim buf3()()() As Double ' filter input, k frames
Dim buf4()() As Double ' filter output, fft input
Dim buf5()() As Double ' fft output
Dim taps()()() As Double ' copies of the filter taps
Dim fftwhin(), fftwhout() As GCHandle
Dim fftwplan() As IntPtr
Dim info() As frameinfo
Public pollinterval As Integer = 20
Public WindowFlag As Boolean = False
''' Window function
Dim win As WindowFunctions
Public Sub New(ByVal Nin As Integer, ByVal threadcountin As Integer, _
ByVal inttime As Double, Optional ByVal kin As Integer = 1, _
Optional ByVal bps As Integer = 4, _
Optional ByVal tapspath As String = Nothing)
N = Nin
threadcount = threadcountin
k = kin
bytespersample = bps
integrationtime = inttime
' loading the filter taps
If k > 1 Then
If N > 2048 Then
Throw (New ArgumentException(sprintf("%s.New - Fourier frame %d is too large.", Me.GetType.Name, N)))
Else
Dim filename As String = sprintf("LPF%06d%02d.fcf", N, k)
If tapspath IsNot Nothing Then
filename = My.Computer.FileSystem.CombinePath(tapspath, filename)
End If
filtertaps = LoadFilterTaps(filename, k, N)
End If
End If
' allocate space for the threads
' create arrays of arrays
ReDim buf2(threadcount - 1)
ReDim buf3(threadcount - 1)
For i As Integer = 0 To threadcount - 1
ReDim buf3(i)(k - 1)
Next
ReDim buf4(threadcount - 1)
ReDim buf5(threadcount - 1)
ReDim taps(threadcount - 1)
For i As Integer = 0 To threadcount - 1
ReDim taps(i)(k - 1)
Next
ReDim fftwhin(threadcount - 1)
ReDim fftwhout(threadcount - 1)
ReDim fftwplan(threadcount - 1)
ReDim info(threadcount - 1) ' value type
' create new arrays for each thread
Dim N2 As Integer = 2 * N
For i As Integer = 0 To threadcount - 1
ReDim buf2(i)(N2 - 1)
For j As Integer = 0 To k - 1
ReDim buf3(i)(j)(N2 - 1)
Next
ReDim buf4(i)(N2 - 1)
ReDim buf5(i)(N2 - 1)
' copy of filter taps for thread i
If k > 1 Then
For j As Integer = 0 To k - 1
ReDim taps(i)(j)(N - 1)
Array.Copy(filtertaps(j), taps(i)(j), N)
Next
End If
' FFTW: handles to the pinned in and out arrays
fftwhin(i) = GCHandle.Alloc(buf4(i), GCHandleType.Pinned)
fftwhout(i) = GCHandle.Alloc(buf5(i), GCHandleType.Pinned)
fftwplan(i) = fftw.dft_1d(N, _
fftwhin(i).AddrOfPinnedObject(), fftwhout(i).AddrOfPinnedObject(), _
fftw_direction.Forward, fftw_flags.Measure)
Next
' window functions from Wikipedia, used when PFB is not
win = New WindowFunctions(N, WindowFunctions.wfType.Blackman)
' number of processor cores
corecount = 0
For Each item As ManagementBaseObject In New System.Management.ManagementObjectSearcher("Select * from Win32_Processor").Get()
corecount += Val(item("NumberOfCores").ToString())
Next
End Sub
Function LoadFilterTaps(ByVal filename As String, ByVal k As Integer, ByVal N As Integer) As Double()()
'read the filter file
Dim filetext As String = My.Computer.FileSystem.ReadAllText(filename, System.Text.Encoding.ASCII)
Dim filelines() As String = filetext.Split(New Char() {Chr(10)})
If filelines.Length - 17 < 4 * N Then
Throw (New ArgumentException(sprintf( _
"%s.process - incorrect number %d of taps found in the filter file %s; %d expected.", _
Me.GetType.Name, filelines.Length - 17, 4 * N)) _
)
End If
' convert the lines of the file
Dim filter(k - 1)() As Double
Dim linecount As Integer = 0
' advance to the first tap in the file
Do Until filelines(linecount).StartsWith("Numerator:") : linecount += 1 : Loop
Dim m As Double = -1
For i As Integer = 0 To k - 1
ReDim filter(i)(N - 1)
Dim f() As Double = filter(i)
For j As Integer = 0 To N - 1
' convert the next filter tap
linecount += 1
f(j) = CDbl(filelines(linecount))
' find the max abs value of the filter taps
If Math.Abs(f(j)) > m Then m = f(j)
Next
Next
' divide by the max
For i As Integer = 0 To k - 1
Dim f() As Double = filter(i)
For j As Integer = 0 To N - 1
f(j) /= m / 1.5
Next
Next
Return filter
End Function
Private Sub process(ByVal info As Object)
' time the computation
Dim s As DateTime = Now
Dim Nminus1 As Integer = N - 1
Dim N2 As Integer = N << 1
Dim N2minus1 As Integer = N2 - 1
Dim sqrtN As Double = Sqrt(CDbl(N))
Dim blocks As Integer
Dim currentframe As Integer = 0
Dim bytesperframe As Integer = bytespersample * N
Dim buf2() As Short = info.buf2
Dim buf4() As Double = info.buf4
Dim buf5() As Double = info.buf5
Dim buf6() As Double = info.buf6
Array.Clear(buf6, 0, buf6.Length)
' preload the filter pipeline
If info.buf1a IsNot Nothing Then
Dim buf1a() As Byte = info.buf1a
blocks = buf1a.Length \ (bytespersample * N)
For i As Integer = 0 To k - 2
' copy a block to a short array
BlockCopy(buf1a, currentframe * bytesperframe, info.buf2, 0, bytesperframe)
' convert block from integer to real
Dim buf3() As Double = info.buf3(currentframe)
For j As Integer = 0 To N2 - 1
buf3(j) = CDbl(buf2(j))
Next
currentframe += 1 ' next frame
Next
Else
' take data from buf1b
Dim buf1b() As Double = info.buf1b
blocks = buf1b.Length \ N2
For i As Integer = 0 To k - 2
Array.Copy(buf1b, currentframe * N2, info.buf3(currentframe), 0, N2)
currentframe += 1 ' next frame
Next
End If
' loop over the blocks
Dim blockcount As Integer = 0
Dim buf1aptr As Integer = currentframe * bytesperframe
For i As Integer = 0 To blocks - k + 1 - 1
Dim bufloc() As Double
If k > 1 Then bufloc = info.buf3(currentframe) Else bufloc = buf4
If info.buf1a IsNot Nothing Then
' copy a block to a short array
BlockCopy(info.buf1a, buf1aptr, buf2, 0, bytesperframe)
For j As Integer = 0 To N2minus1
bufloc(j) = CDbl(buf2(j))
Next
Else
' take data from buf1b
'Array.Copy(info.buf1b, currentframe * N2, bufloc, 0, N2)
Dim N16 As Integer = N << 4
BlockCopy(info.buf1b, i * N16, bufloc, 0, N16)
End If
If k > 1 Then
' filter the block
Array.Clear(buf4, 0, N2)
Dim buf3l(), taps() As Double
For j As Integer = 0 To k - 1
Dim cfj As Integer = (currentframe + j + 1) Mod k
Dim l2 As Integer = 0
buf3l = info.buf3(cfj)
taps = info.taps(j)
For l As Integer = 0 To Nminus1
Dim tap As Double = taps(l)
' real part
buf4(l2) += tap * buf3l(l2) : l2 += 1
' imaginary part
buf4(l2) += tap * buf3l(l2) : l2 += 1
Next
Next
Else
' apply window function if enabled
If WindowFlag Then win.mult2(bufloc)
End If
' increment the frame number
currentframe = (currentframe + 1) Mod k
' compute the fourier transform: buf4 -> buf5
fftwf.execute(info.fftwplan)
' accumulate the spectral density
For j As Integer = 0 To Nminus1
Dim idx As Integer = j << 1
Dim a As Double = buf5(idx)
Dim b As Double = buf5(idx + 1)
buf6(j) += a * a + b * b ' faster
Next
' quit if the time limit is exceeded
If i Mod pollinterval = 0 AndAlso Now.Subtract(s).TotalSeconds > info.timelimit Then
blockcount = i + 1
Exit For
End If
Next
' divide by the number of frames processed in the time limit
If blockcount = 0 Then blockcount = blocks - k + 1
Dim scale As Double = 1 / (blockcount * N)
For i As Integer = 0 To nminus1
buf6(i) *= scale
Next
' print performance data
'Static count As Integer = 0
'Dim interval As Double = Now.Subtract(s).TotalMilliseconds
'count = count + 1
'printf("time %6.2f ms, %5.1f us/block, blocks processed %d/%d (%4.1f percent)", _
' interval, interval / blockcount * 1000.0, blockcount, blocks, blockcount / blocks * 100)
End Sub
'''
'''
'''
''' Raw byte data to be interpreted as interleaved I and Q shorts.
'''
''' Arrays for receiver data and output spectra are provided by the caller.
Public Function run(ByVal data() As Byte, ByVal result() As Double) As Double()
' info structure going to the thread
info(currentthread).buf1a = data ' which receiver buffer to use
info(currentthread).buf1b = Nothing ' no double data
' intra-class assignments
info(currentthread).buf2 = buf2(currentthread)
info(currentthread).buf3 = buf3(currentthread) ' fft input
If k > 1 Then
info(currentthread).taps = taps(currentthread) ' filter taps
Else
info(currentthread).taps = Nothing ' not used
End If
info(currentthread).buf4 = buf4(currentthread) ' fft output
info(currentthread).buf5 = buf5(currentthread)
info(currentthread).buf6 = result ' spectrum returned
info(currentthread).fftwplan = fftwplan(currentthread) ' fftw plan to be run by the thread
info(currentthread).timelimit = integrationtime * (Min(corecount, threadcount) * 0.85 - 0.4)
' start the thread
ThreadPool.QueueUserWorkItem(AddressOf process, info(currentthread))
' increment the current thread
currentthread = (currentthread + 1) Mod threadcount
' return the spectrum of the oldest thread
Return info(currentthread).buf6
End Function
Public Function run(ByVal data() As Double, ByVal result() As Double) As Double()
' structure going to the thread
info(currentthread).buf1a = Nothing ' no byte data
info(currentthread).buf1b = data ' double data
' intra-class assignments
info(currentthread).buf2 = Nothing ' not used
info(currentthread).buf3 = buf3(currentthread) ' fft input
' filter taps
If k > 1 Then
info(currentthread).taps = taps(currentthread) ' filter taps
Else
info(currentthread).taps = Nothing ' not used
End If
info(currentthread).buf4 = buf4(currentthread) ' fft output
info(currentthread).buf5 = buf5(currentthread)
info(currentthread).buf6 = result ' spectrum returned
info(currentthread).fftwplan = fftwplan(currentthread) ' fftw plan to be run by the thread
info(currentthread).timelimit = integrationtime * (Min(corecount, threadcount) * 0.85 - 0.4)
' queue the thread
ThreadPool.QueueUserWorkItem(AddressOf process, info(currentthread))
' increment the current thread
currentthread = (currentthread + 1) Mod threadcount
' return the spectrum of the oldest thread
Return info(currentthread).buf6
End Function
''' Stucture used for providing information to threads used for computing.
''' Provide data in buf1a or buf1b, but not both. Leave a null pointer in the other.
Public Structure frameinfo
Dim corecount As Integer
Dim buf1a() As Byte ' alternating I and Q shorts;
Dim buf1b() As Double ' alternating I and Q doubles
Dim buf2() As Short
Dim buf3()() As Double
Dim buf4() As Double
Dim buf5() As Double
Dim buf6() As Double ' spectrum
Dim taps()() As Double ' filter taps
Dim fftwplan As IntPtr ' FFTW plan
Dim timelimit As Double
End Structure
Shadows Sub finalize()
close()
MyBase.Finalize()
End Sub
Public Sub close()
' cleanup
waittofinish()
For i As Integer = 0 To threadcount - 1
fftw.destroy_plan(fftwplan(i))
fftwhin(i).Free()
fftwhout(i).Free()
Next
End Sub
Public ReadOnly Property Latency()
Get
Return threadcount
End Get
End Property
Public Sub waittofinish()
Sleep((threadcount * (integrationtime + 0.03)) * 1000)
End Sub
Delegate Sub elog(ByVal a As String, ByVal b As String, ByVal c As Boolean, ByVal d As System.Text.Encoding)
''' Write exception information to a log file.
''' The exception to be logged.
Public Shared Sub log(ByVal ex As Exception, Optional ByVal filename As String = "SpectralProcessing.log")
Static e As elog = AddressOf My.Computer.FileSystem.WriteAllText
Static a As System.Text.Encoding = System.Text.Encoding.ASCII
e(filename, "=========================================" & vbNewLine, True, a)
e(filename, ex.Message & vbNewLine, True, a)
e(filename, ex.StackTrace & vbNewLine, True, a)
If ex.InnerException IsNot Nothing Then
e(filename, "= = = = = = = = = = = = = = = = = = = = =" & vbNewLine, True, a)
e(filename, ex.InnerException.StackTrace & vbNewLine, True, a)
End If
End Sub
End Class