I missed one crucial fact when I looked at f2py and that is that it requires the source as it compiles that source as part of generating the py wrapper.
It still seems possible to call from python though writing a wrapper from scratch will be a little more involved (pyfort seems to be one option, doing it manually is the other)
It still seems possible to call from python though writing a wrapper from scratch will be a little more involved (pyfort seems to be one option, doing it manually is the other)
I'm just curious does loading the .dll and calling a function with the proper attention to the types and byref, etc. not work with Fortran? I thought Numpy had all the slice and dice array stuff to go between the two built-in. Sorry if I'm being ignorant, I'm sure you know better than me.
Well I've never called a fortran dll from python before so this is new to me too
My interpretation of Using Python as glue — NumPy v1.14 Manual is
- if you have the source, you can use f2py to ultimately create a "friendly" interface with good performance characteristics in an extension module
- if you don't have the source then an creating an extension module is probably not possible due to the init requirements (How to extend NumPy — NumPy v1.14 Manual)
- ctypes is the mechanism for accessing arbitrary DLLs and numpy does provide a mechanism for converting an ndarray to a ctype type so it seems this should ultimately work once you arrive at the correct incantation
The apparent requirement to use ctypes is why I said it seems like it will be a little more involved.
That's how we do the interface to our USB module, but our dll is in C. In the system module there is a function to load the dll and make all the function names visible (mydll.myfunction())
The problem is how the complex arrays are passed in Fortran. My experience is very old but when translating Fortran to C for an FFT library IIRC the complex type was two separate N x 1 arrays (for 1 dimension), Python's built in complex type is a true single object but I don't see what could be easier than passing real() and imag() parts separately.
This is all just bookkeeping in my view. The worst case I ever dealt with was packed 24bit audio (odd bytes, 3) going from a little endian to big endian processor where the bytes are swapped on 2 byte boundaries.
It's just a boring thing to write and annoying to debug
Yes it was
I was just looking at the no source avaible aspect of this problem.Simple it is.
Yes, with a fully automated system it is OK to have whatever is mechanically simple and just let the CNC work out the correct combination of co-ordinate moves.
For a manually scanned system it is more important to make it simple to move the system to the next measurement.
It would be convenient if the "North pole" was on the speaker axis.
The mechanics I proposed as my first idea would require the speaker to be on it's back, pointed up.
That's not ideal, it puts the important measurement axis on the usually small dimension of the room.
Also an inconvenience for speakers like the "Voice of the Theater" A7 style clones I used to own, with the horns just placed on the top.
This style of speaker would also be inconvenient for Earl's idea to place them on their sides to measure verticals.
One option would be to rotate my idea so the North Pole axis (theta axis) is horizontal.
That would require a strut out the back of the speaker platform but would have the benefit that we wouldn't have to put the substantial speaker mount shaft thru the middle of the rotation system.
Instead we could use a standard stepper motor and gearhead, counter balanced, the balance could even be a mass damper to steady the movement of the microphone after it steps to it's new position.
I did see that one, but when you said a "paper" I took it to mean a paper in a journal or website rather than an actual piece of paper.
Not sure what that says about me or my expectations.
Anyway, now it's clarified and I have reread the post.
Is the assumed source radius parameter a a variable that can be input by the user or is a fixed value assumed in the code?
Best wishes
David.
Last edited:
It's a variable in my code. I usually set it to a value that yields the same volume of sphere as the volume of the box, but in reality what it is set at is not significant unless you are doing a reconstruction of the source. That's because it cancels for far field reconstruction as I describe in the "paper" (or whatever you want to call it.)
I think a python wrapper for the FFT call is
Code:
__all__ = ['fft']
import numpy as np
import os
from ctypes import *
# assume the ddl is in the same path
_path = os.path.dirname('__file__')
lib = np.ctypeslib.load_library('MeasCalcs', _path)
func = getattr(lib, 'FFT')
func.restype = None
func.argtypes = [np.ctypeslib.ndpointer(float, flags='aligned, contiguous'), POINTER(c_int32)]
# assumes input is an 1 dim ndarray (i.e. should validate this before the func call)
def fft(data):
data = np.require(np.asanyarray(data, float, ['CONTIGUOUS', 'ALIGNED'])
len = pointer(c_int32(data.shape[0]))
func(data, len)
return data.contentsI am assuming the FFT subroutine changes the contents of the data_ array (though there may be some additional flags required to make that happen, tbd once I test it on the actual function)
Last edited:
It is confusing I hope Earl can help explain. Does it maybe work like np.fft.rfft?
I described this before somewhere.
The FFT takes in a real array of N and passes out a complex array of N/2 in the same array space that was passed in. This WAS commonplace in FORTRAN from the days when memory was at a premium. Today one would just pass in one array and pass out another.
As I said before, the FFT is an array of N/2 + 1, and the zeroth element is the real part of the first element in the array passed out and the real part of the N/2+1 element is is the imaginary part of the zeroth element. This is possible because the imaginary parts of both these elements has to be real. Again, commonplace in the "OLD" days.
These are precisely the kinds of complications that will make blending very old code with a modern language. Things are just done differently now.
The FFT takes in a real array of N and passes out a complex array of N/2 in the same array space that was passed in. This WAS commonplace in FORTRAN from the days when memory was at a premium. Today one would just pass in one array and pass out another.
As I said before, the FFT is an array of N/2 + 1, and the zeroth element is the real part of the first element in the array passed out and the real part of the N/2+1 element is is the imaginary part of the zeroth element. This is possible because the imaginary parts of both these elements has to be real. Again, commonplace in the "OLD" days.
These are precisely the kinds of complications that will make blending very old code with a modern language. Things are just done differently now.
Last edited:
Not that differently. I think just a tiny bit of bookkeeping difference from the built it real fft function. You are right no one thinks twice these days about making multiple copies of even millions of element floating point arrays.
right I see, numpy supports this via numpy.ndarray.view — NumPy v1.14 Manual
i.e. I can just return output.view(dtype=np.complex128) to return an array of complex numbers instead
FWIW it was the fact it was described as intent(in) rather than intent(inout) that made me think something might be missing.
i.e. I can just return output.view(dtype=np.complex128) to return an array of complex numbers instead
FWIW it was the fact it was described as intent(in) rather than intent(inout) that made me think something might be missing.
Whatever that code was that you showed may be powerful, but it is very cryptic, not at all easy to see what it means.
Yea, weak programming on my part. "Intent" should not have been used.
Sorry Earl I won't bother in the future, I simply showed an 8 element real array has a 5 element complex (N/2 +1) FFT and goes back again to a real array with the inverse transform. No variable typing necessary at all.
I'll leave you guys to do this.
You should probably leave it here since that is one purpose of the thread and there are not any other takers. I would gladly help in any way I can as limited as that may be.