Klippel Near Field Scanner on a Shoestring

[3ll3d00d]

See, I can't follow that as there is just not enough information.

take it from me that it is the same, for reference numpy.logspace — NumPy v1.14 Manual

I just take the data and whereever the window is set to place data and beyond that zeros. I don't see what a power of two has to do with it.

your fft function accepts power of 2 sized arrays only. Does your code always pass in the same sized array (and if so, what size is it?) or do just you use the next biggest power of 2?

I don't know "tukey" as a window

tukey is a named window function provided by scipy, there is a description along with references and images showing the window in scipy.signal.windows.tukey — SciPy v1.1.0 Reference Guide

this is described as a tapered cosine window and "can be regarded as a cosine lobe of width αN/2 that is convolved with a rectangular window of width (1 − α/2)N"

I can set α though the default is 0.5

if this isn't suitable then I can implement something else if you (or anyone else) can specify it

I can't follow your explanations about the arrays. I don't understand

A[2:].view(complex128) == X[2:257]

but to me that is 257 points which is correct

.view(complex128) basically just changes the datatype presented by the array without changing the underlying data

for example given a 4 element float64 array [1.0,0.0,1.0,0.0] then view(complex128) will return a 2 element complex128 array [1.0+0.0i,1.0+0.0i]

i.e. each 16 consecutive bytes (aka 2 array slots) in the input becomes a single complex number in the output

 

[3ll3d00d]

@gedlee is this linToLog problem a blocker for continuing further or can we develop the next steps of the UI in parallel working on the assumption that we will solve the linToLog problem at some point?

the reasons to proceed are twofold

1) bashing away at dll problems is extremely tedious and I don't have any real clues as to what is wrong atm
2) making use of some of the other functions in the dll may shed some light on what is wrong with linToLog

however if the next steps simply don't work without that linToLog output then obviously we have to pause here and work this out

 

[gedlee]

Everything after LogToLin assumes the data is linearly spaced. I highly doubt that we could continue without this step.

For now we can set aside the window and just null out all data past the window line out to the power of 2 number of samples assumed in the code.

But here is how I do the window:

    Structure Window
        Dim On_ As Boolean
        Dim Start As Single  ' in ms.
        Dim StartSet As Boolean
        Dim End_ As Single   ' in ms.
        Dim EndSet As Boolean
        Function intTimeStart() As Integer  '  integer value in data set
            intTimeStart = CInt(Start * SampleRate_ / MilliSec)
        End Function
        Function intTimeEnd() As Integer '  integer value in data set
            intTimeEnd = CInt(End_ * SampleRate_ / MilliSec)
        End Function
    End Structure

    Public win As Window

   Function WindowData(ByVal current As Integer) As Double
        Dim loc As Integer
        WindowData = 1.0
        If current < (win.intTimeEnd - BlendWidth_ * sampleRateMS) Then
            WindowData = 1.0
        ElseIf current > (win.intTimeEnd + BlendWidth_ * sampleRateMS) Then ' above blend
            WindowData = 0.0
        Else        ' blend
            loc = CInt(current - win.intTimeEnd + BlendWidth_ * sampleRateMS)
            WindowData = (Math.Cos(loc * Math.PI / (BlendWidth_ * sampleRateMS) / 2.0) + 1.0) / 2.0
        End If

    End Function

"Blendwidth" is a variable that is usually about 1 ms.

Like I said, debugging DLLs has been the most difficult thing for me in the past and I cannot see how one debugs them if you cannot debug through the call, as I can do in Visual Studio.

 

[3ll3d00d]

Like I said, debugging DLLs has been the most difficult thing for me in the past and I cannot see how one debugs them if you cannot debug through the call, as I can do in Visual Studio.

it would certainly be useful in this case

anyway I thought that I would attempt to rule out python so I created a minimal, complete and verifiable example using C# in VS. I ran this with .net framework 4.6.1 and also 4.5.

The input.txt file mentioned was generated by the python app to check the same data is used in each case and is attached.

View attachment input.txt

I tried this in both VS2013 and VS2017 with same result by adding the path to the required DLLs to the PATH, creating a new consoleapp project, adding a single additional reference to System.Numerics and switching the target to x86.

Hopefully this is sufficient for further debugging or I can share the vs project instead.

using System;
using System.Runtime.InteropServices;
using System.IO;
using System.Numerics;

namespace ConsoleApp2
{
    class Program
    {
        static void Main(string[] args)
        {
            // read the windowed input data into a double array
            var lines = File.ReadAllLines("C:\\Users\\mattk\\source\\repos\\ConsoleApp3\\ConsoleApp3\\input.txt");
            var n = lines.Length;
            var vals = new double[n];
            for (var i = 0; i < n; i++)
            {
                vals[i] = Double.Parse(lines[i]);
            }
            Console.WriteLine("Parsed " + n + " samples");
            // FFT it
            MeasCalcs.FFT(vals, ref n);
            // unwrap the output into n/2+1 complex values
            var complexN = (n / 2) + 1;
            var fftData = new Complex[complexN];
            fftData[0] = new Complex(vals[0], 0);
            fftData[complexN - 1] = new Complex(vals[1], 0);
            for (var i = 2; i < n; i += 2)
            {
                fftData[i / 2] = new Complex(vals[i], vals[i + 1]);
            }
            for (var i = 0; i < fftData.Length; i++)
            {
                Console.WriteLine("fft " + i + ":" + fftData[i]);
            }
            // create a log spaced array of frequencies from 20-20k
            var logSteps = 128;
            var logFreqs = new double[logSteps];
            for (int i = 0; i < logSteps; i++)
            {
                logFreqs[i] = 2.0 * Math.Pow(10.0, 1.0 + (3.0 * (i / ((double)logSteps - 1.0))));
                Console.WriteLine("logFreq: " + i + " - " + logFreqs[i]);
            }
            // initialise an output array with zeroes
            var linToLogData = new Complex[logSteps];
            for (int i = 0; i < logSteps; i++)
            {
                linToLogData[i] = new Complex(0, 0);
            }
            Console.WriteLine(">> LinToLog");
            double freqStep = 48000 / n;
            int fftPts = fftData.Length;
            MeasCalcs.LinToLog(fftData, freqStep, fftPts, linToLogData, logSteps, logFreqs);
            Console.WriteLine("<< LinToLog");
            // dump as an FRD
            Console.WriteLine("freq mag phase");
            for (var i = 0; i < linToLogData.Length; i++)
            {
                Console.WriteLine(logFreqs[i] + " " + linToLogData[i].Magnitude + " " + linToLogData[i].Phase);
            }
            Console.WriteLine(" Done!");
        }
    }

    class MeasCalcs
    {
        [DllImport("MeasCalcs.dll", CallingConvention = CallingConvention.Cdecl)]
        public static extern void FFT(double[] vals, ref int size);

        [DllImport("MeasCalcs.dll", CallingConvention = CallingConvention.Cdecl)]
        public static extern void LinToLog(Complex[] dataIn, double deltaF, int fftPts, Complex[] dataOut, int logPts, double[] logFreqs);
    }
}

 

[gedlee]

First congrats to England on making the quarters. My bet is now on Belgium though.

I can probably run your example above, except that I am missing the referenced "system.numerics". Where is that from? It's not in my .net 4.6. My VS does not recognize it. With it, I should be able to debug through the call and get some answers. I am completely comfortable with this C# code. If we can work through problems in this language and then you can translate into Python if you want and we should be OK. But Python is just too foreign for me and I do not want to learn another all-new language - sorry.

 

[3ll3d00d]

belgium and france definitely look dangerous, only needs 2 wins now though!

re numerics, I found it as per c# - How to add a reference to System.Numerics.dll - Stack Overflow

on my installation it lives in C:\Program Files (x86)\Reference Assemblies\Microsoft\Framework\.NETFramework\v4.5\System.Numerics.dll

I believe this was added in .NET 4.0 so you should have it. There are a number of Qs like "where is the numerics dll" on SO though, not sure why that happens. Try the instructions in the link above anyway and see if that helps.

no problem with writing examples in c# as necessary to work this out (I just choose not to write in a language like c# in my spare time)

 

[gedlee]

Yea, that was easy. Don't understand why it was not the default and had to be added separately.

 

[gedlee]

Ok, I could compile and run the code, but I am currently having what I call "permissions hell". Somehow I have my "permissions" all screwed up and can't get it sorted out. I use 8.1 and I may have to upgrade to 10, which will take me quite awhile to get back to productive. The code tells me that it can't open "input.txt" because I don't have permission, but yet it's in my own "\users|" directory and the file shows that "everyone" have "full control". I just don;t get "permissions" in the windows environment. When it gets screwed up things go really bad.

No other computer on my network can access the supposed shares on my computer and I can't go the other way either. Its as if all communication between my computer and any other is locked out and now even I am locked out of certain places. I don't get it.

 

[gedlee]

I was able to debug through the call, but am getting an error that I have never seen before. The code seems to run fine within the FORTRAN code, but fails on return. I get an error in the DLL that says

Unhandled exception at 0x77f76214 in DLL_call_test.exe: 0xC0000374: A heap has been corrupted.
The Common Language Runtime cannot stop at this exception. Common causes include: incorrect COM interop marshalling, and memory corruption. To investigate further, using native-only debugging.

and

		YPP	CXX0072: Error: type information missing or unknown
		FREQIN	CXX0072: Error: type information missing or unknown

Both of which are declared. FreqIn is passed in and likely the culprit as it can't pass out of the "type" is unknown. This may be a C# thing as this is a new error to me.

But this has been a big jump! That's for sure. Being able to debug through the call is essential to debugging inter-language calls.

 

[3ll3d00d]

interesting

is FreqIn passed in?

the doc you posted earlier said

subroutine LinToLog(    DataIn,     &   ! FFT data in at all angles  
                            DeltaF,     &   ! the frequency delta
                            NumFFTPts,  &   ! the number of FFT points
                            Dataout,    &   ! The interpolated log frequency output
                            NumLogPts,  &   ! the number of data points being output
                            Freqs       )   ! the Array of output frequencies
                         
     
!DEC$ ATTRIBUTES REFERENCE :: DataIn
!DEC$ ATTRIBUTES Value :: NumFFTPts, DeltaF
!DEC$ ATTRIBUTES REFERENCE :: Dataout
!DEC$ ATTRIBUTES REFERENCE :: Freqs
!DEC$ ATTRIBUTES Value :: NumLogPts
integer, intent(in):: NumFFTPts                             ! input array size constants
integer, intent(in) :: NumLogPts                                        ! number of
smoothed points ‐ log scaled
COMPLEX(8), intent(in):: DataIn(0:NumFFTPts)             ! input frequency data  !
changed
COMPLEX(8), intent(out) :: DataOut(0:NumLogPts)         ! output data array ! changed
real(8), intent(in) :: DeltaF
real(8), intent(in) :: Freqs(0:NumLogPts)
real(8)  Freqin(0:NumFFTPts)

I took this to mean Freqin was a local variable and you're populating this based on the DeltaF value of (i.e. it will contain i*DeltaF)

if so, does that suggest the value for DeltaF is bad perhaps?

 

[mkc]

Just a thought as this bumped my memory. Perhaps not the issue here and I haven't looked into this in detail.

However, I seem to remember in a project (more than 10 years ago), that using an DLL (in this case written in c++) from C#, that the C# garbage collector could be a little bit too effective if the memory was passed on to un-managed code. It would just de-allocate the memory.

Again, it was a long time ago so my memory of the details are not great. Just wanted to throw it in here to say that this could result in strange behavior.

Cheers,
Mogens

 

[3ll3d00d]

It's a good point but I don't think this is an issue for 2 reasons

1) the memory is allocated on the c# (or python) side and then passed to the dll, this means the objects are not eligible for collection
2) exact same behaviour is seen in 2 completely different runtimes (python and c#)

the other oddity is that it doesn't happen every single time, it does happen most of the time but it still intermittently succeeds albeit mostly after I've left the computer for a while. This could be coincidence or it could be a significant indicator towards what is going on. Hindsight will probably tell us the answer to that one.

finally I'm not 100% sure how that Complex object is passed to a dll, it does say it is a struct not a class so I assume this will work as expected but I can't say I'm 100% sure about that.

 

[gedlee]

One problem that I found is that NumLogPts and NumFFTpts should really be MaxLogPt and MaxFFTpt in your code as inside the FORTRAN all the arrays are one too large (ie NumLogPts = MaxLogPt + 1 given the zero base). But even correcting this in your code, I get an internal error on a call to a routine that I did not write. Its from Numerical Recipes, "splint". I did modify it for complex numbers but it has still always worked in my program. The value of "NumFFTPts" becomes undefined in the subroutine splint and then it blows up. This may also be the index being off by one, I need to look more closely at it. Hence, at this point there is nothing you can do except maybe relook at you code in light of the error that I described above (my bad!) I'm sure that this all has to do with getting the sizes of the arrays correct.

 

[gedlee]

I also noticed that in the FORTRAN code Freqs(128) is NaN

 

[3ll3d00d]

a quick test with simply reducing those values by 1 results in something that runs consistently (albeit see previous post for intermittent successes)

let me try the python side....

 

[3ll3d00d]

seems to work ok in python as well now, thanks for figuring that one out, so I guess we can move on from here now?

one question before that though concern the amplitude of the output from the FFT routine... IIRC if you want the amplitude out of the fft to equal the input then you need to do something like normalise based on the no of pts in the input (I am thinking of something like amplitude of (signal) after FFT operation? -
MATLAB Answers - MATLAB Central). Does your FFT routine handle this or is something like this something I need to do?

secondly, for plotting the FR, is there a preferred way to scale the amplitude? relative dB (where 1 = 0dB assuming a peak value of 1 in the input data)? something else?

beyond that, what are the next steps in building this out? I have some UI bits and pieces I could polish but prefer to hack out something functional first before doing any polishing.

 

[gedlee]

This will be a problem.

I have never had a calibrated system, which requires calibration from mic all the way through to the results. So I have never bothered with calibration in the code for things like the FFT. Absolute values of SPL are kind of unimportant to me, so I just normalize the results at the end to be 0 dB average over about 1/3 octave at 1 kHz. I suggest we not worry about calibration at this point at later go through and make that happen.

I see things this way. Let's get something simple and usable and then see if we can't get some funding through the various means to make a really competitive system. I am not interested in any compensation, but the potential for benefit to others might make things go more smoothly. I don't think that we should ever sell the software, but we could ask for donations from users. If people will pay Klippel's prices we should be able to make something competitive worth some donations.

 

[Kravchenko_Audio]

You can't work with a basic microphone sensitivity in millivolt/Pascal?

It's at the least something that is a baseline that is easily found in almost all microphones.

 

[scott wurcer]

This might be a little off topic but of interest to those that might want to incorporate the acquisition of the measurement data into a single Python program. SoX is very popular for audio file manipulation, recording, and playing. I know under Windows Python has a pyaudio interface to the free PortAudio implementation but some searching found an already built-in interface (subprocess) that works with SoX making all its capabilities available. Sox (SoX - Sound eXchange | HomePage) is a one click install for those that want it easy. If you have Anaconda Python and SoX with a Windows path entry to sox.exe there is nothing else to do.

from scipy.io.wavfile import read, write
import matplotlib.pyplot as plt
import numpy as np
import subprocess
import time

#edit here to add a HOME directory, file paths, etc.
RECFILE = 'C:\Users\scott\Desktop\\tempau2.wav'
PLAYFILE = 'C:\Users\scott\Desktop\\tone2.wav'
#end edit

# Build the argument strings for SoX. The commands must self terminate, play ends at
# the end of file and record ends by using the trim option, trim 0 4 simply records
# 0 to 4 seconds and stops. You can't cont-C to stop recording here.

args_play = ['sox','-q',PLAYFILE,'-t','waveaudio','default']
args_rec = ['sox','-q','-t','waveaudio','default','-t','wavpcm',RECFILE,'trim','0','4']

# Force SoX to run without opening a cmd window, the -q option in the args
# also forces SoX into quiet mode that minimizes latency

startupinfo = subprocess.STARTUPINFO()
startupinfo.dwFlags |= subprocess.STARTF_USESHOWWINDOW

# Simultaneous play and record from seperate processes usually works these days.
# My Scarlett 2i2 didn't (after hours of aggravation I simply checked and it does not
# with any application (sans ASIO) under Win 10). A cheap Soundblaster USB worked fine. You will need to set
# default devices in the Windows control panel and fiddle with the mixer settings
# to get levels right, 48k at 16 bits worked fine. Don't be fooled into thinking you are
# getting 24 bits under Windows, none of the free PortAudio implementations (AFAIK) support it, in
# general ASIO is required for this. Trust me under Win (not Linux) you can really be fooled.

sox_r = subprocess.Popen(args_rec, startupinfo = startupinfo)
sox_p = subprocess.Popen(args_play, startupinfo = startupinfo)

# Python keeps going, I found simply sleeping until the played file is over works since
# I could not find a foolproof process to signal completion under Windows. The '5' below
# is the length of the played file, a general solution would test for the length and
# insert it. The recording needs to end first with enough time to close the file, otherwise
# there can be garbage in the file. Padding the files and accounting for this should not
# be a problem.

time.sleep(5)

The picture is a loopback of a 1kHz tone at around 30,000 16 bit codes.

 

[scott wurcer]

You can't work with a basic microphone sensitivity in millivolt/Pascal?

It's at the least something that is a baseline that is easily found in almost all microphones.

You should be able to run a simple cal cycle using the above with the mic input simply knowing the sensitivity of the mic. This is something that could be added but most useful information is based on relative data, there are not many cases where absolute SPL plus minus several dB really matters.