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