arikitik
Янв. 6, 2010 00:07:12
Есть файл, строгой структуры, в котором записан, в том числе массив данных в формате 4-byte IBM floating point. При помощи модуля struct я получаю значения, но они отличаются от действительных. При визуализации получается график такого же вида, как правильный, но амплитуды другие. Подскажите, в чём может быть проблема.
фрагмент кода
fname = sys.argv
sgy = open(fname, ‘rb’)
sgy.seek(3212)
traces = struct.unpack('>h', sgy.read(2))
sgy.seek(3220)
samples = struct.unpack('>h',sgy.read(2))
data = np.zeros(samples*traces)
data = data.reshape(120,5000)
sgy.seek(3600)
for x in range(traces):
sgy.seek(sgy.tell()+240)
for y in range(samples):
data = struct.unpack('>f',sgy.read(4))
sgy.close()
Андрей Светлов
Янв. 6, 2010 20:33:36
Насколько я смог понять, IBM floating point несколько отличается от IEEE-754 floating point (последняя используется сопроцессором вообще и Питоном в частности).
Для 32 битного числа имеем:
IEEE-754: Sign bit: 1 bit, Exponent width: 8 bits, Significand precision: 24 (23 explicitly stored)
IBM: Sign - 1bit, exponent: 6 bits, mantissa: 25 bit
Экспонента на 2 бита короче, а мантиса соответственно длиннее.
Модуль struct вам не поможет - ищите или пишите преобразование специально для IBM формата.
arikitik
Янв. 7, 2010 02:08:18
печально(( спасибо за разъяснение!
Андрей Светлов
Янв. 7, 2010 04:18:04
Не очень-то и печально.
Смотрите:
- ibm формат
http://en.wikipedia.org/wiki/IBM_Floating_Point_Architecture - аналогичный IEEE-754
http://en.wikipedia.org/wiki/Single_precision_floating-point_formatПитон работает с double (процессор, грубо говоря, тоже). Т.е. можно развернуть IBM-32 в IEEE double без потери точности
Это - наш стандарт:
http://en.wikipedia.org/wiki/Double_precision_floating-point_formatберем и побитово переносим из ibm в ieee. Экспонента дополняется старшим битом слева, а мантисса нулем справа.
Когда все заработает - можно быстро перенести код на Cython для более или менее приличной производительности.
Вот кусок на С, который делает то же самое судя по описанию.
http://bytes.com/topic/c/answers/221981-c-code-converting-ibm-370-floating-point-ieee-754-a
arikitik
Янв. 7, 2010 21:12:57
Да, Вы правы) Погуглил на эту тему и всё получилось. Надо бы действительно на Си перенести.
Андрей Светлов
Янв. 7, 2010 22:29:42
Будут проблемы с созданием C extension - спрашивайте. На этом вопросе собаку съел