高效从二进制文件创建Numpy数组的方法

Question

我有一些非常大的数据集，这些数据存储在硬盘上的二进制文件里。下面是文件结构的一个例子：

文件头

149 Byte ASCII Header

记录开始

4 Byte Int - Record Timestamp

样本开始

2 Byte Int - Data Stream 1 Sample
2 Byte Int - Data Stream 2 Sample
2 Byte Int - Data Stream 3 Sample
2 Byte Int - Data Stream 4 Sample

样本结束

每个记录有122,880个样本，每个文件有713个记录。这样总共的大小是700,910,521字节。有时候样本率和记录的数量会有所不同，所以我需要编写代码来检测每个文件中的数量。

目前我用来将这些数据导入数组的代码是这样的：

from time import clock
from numpy import zeros , int16 , int32 , hstack , array , savez
from struct import unpack
from os.path import getsize

start_time = clock()
file_size = getsize(input_file)

with open(input_file,'rb') as openfile:
  input_data = openfile.read()

header = input_data[:149]
record_size = int(header[23:31])
number_of_records = ( file_size - 149 ) / record_size
sample_rate = ( ( record_size - 4 ) / 4 ) / 2

time_series = zeros(0,dtype=int32)
t_series = zeros(0,dtype=int16)
x_series = zeros(0,dtype=int16)
y_series = zeros(0,dtype=int16)
z_series = zeros(0,dtype=int16)

for record in xrange(number_of_records):

  time_stamp = array( unpack( '<l' , input_data[ 149 + (record * record_size) : 149 + (record * record_size) + 4 ] ) , dtype = int32 )
  unpacked_record = unpack( '<' + str(sample_rate * 4) + 'h' , input_data[ 149 + (record * record_size) + 4 : 149 + ( (record + 1) * record_size ) ] ) 

  record_t = zeros(sample_rate , dtype=int16)
  record_x = zeros(sample_rate , dtype=int16)
  record_y = zeros(sample_rate , dtype=int16)
  record_z = zeros(sample_rate , dtype=int16)

  for sample in xrange(sample_rate):

    record_t[sample] = unpacked_record[ ( sample * 4 ) + 0 ]
    record_x[sample] = unpacked_record[ ( sample * 4 ) + 1 ]
    record_y[sample] = unpacked_record[ ( sample * 4 ) + 2 ]
    record_z[sample] = unpacked_record[ ( sample * 4 ) + 3 ]

  time_series = hstack ( ( time_series , time_stamp ) )
  t_series = hstack ( ( t_series , record_t ) )
  x_series = hstack ( ( x_series , record_x ) )
  y_series = hstack ( ( y_series , record_y ) )
  z_series = hstack ( ( z_series , record_z ) )

savez(output_file, t=t_series , x=x_series ,y=y_series, z=z_series, time=time_series)
end_time = clock()
print 'Total Time',end_time - start_time,'seconds'

这个过程目前大约需要250秒来处理一个700MB的文件，对我来说这个时间似乎太长了。有没有更高效的方法可以做到这一点呢？

最终解决方案

使用numpy的fromfile方法，并配合自定义的数据类型，使得运行时间缩短到了9秒，比上面的原始代码快了27倍。最终的代码如下。

from numpy import savez, dtype , fromfile 
from os.path import getsize
from time import clock

start_time = clock()
file_size = getsize(input_file)

openfile = open(input_file,'rb')
header = openfile.read(149)
record_size = int(header[23:31])
number_of_records = ( file_size - 149 ) / record_size
sample_rate = ( ( record_size - 4 ) / 4 ) / 2

record_dtype = dtype( [ ( 'timestamp' , '<i4' ) , ( 'samples' , '<i2' , ( sample_rate , 4 ) ) ] )

data = fromfile(openfile , dtype = record_dtype , count = number_of_records )
time_series = data['timestamp']
t_series = data['samples'][:,:,0].ravel()
x_series = data['samples'][:,:,1].ravel()
y_series = data['samples'][:,:,2].ravel()
z_series = data['samples'][:,:,3].ravel()

savez(output_file, t=t_series , x=x_series ,y=y_series, z=z_series, fid=time_series)

end_time = clock()

print 'It took',end_time - start_time,'seconds'

Answer 1

一个明显的低效问题是，在循环中使用了 hstack：

  time_series = hstack ( ( time_series , time_stamp ) )
  t_series = hstack ( ( t_series , record_t ) )
  x_series = hstack ( ( x_series , record_x ) )
  y_series = hstack ( ( y_series , record_y ) )
  z_series = hstack ( ( z_series , record_z ) )

在每次循环中，这都会为每个系列分配一个稍微大一点的数组，并把到目前为止读取的所有数据复制到这个新数组里。这涉及到很多不必要的数据复制，可能还会导致内存碎片化的问题。

我会把 time_stamp 的值先放到一个列表里，最后再进行一次 hstack，对 record_t 也是一样的做法。

如果这样做还没有明显提升性能，我会把循环里的内容注释掉，然后逐步把每一部分重新加回去，看看具体是在哪个地方花费了时间。

Answer 2

Numpy支持通过 numpy.memmap 将数据直接映射到类似数组的对象中。你可以通过内存映射文件，提取你需要的数据。

为了确保数据的字节顺序正确，只需对你读取的数据使用numpy.byteswap。你可以用一个条件表达式来检查你电脑的字节顺序：

if struct.pack('=f', np.pi) == struct.pack('>f', np.pi):
  # Host is big-endian, in-place conversion
  arrayName.byteswap(True)

Answer 3

一些提示：

• 不要使用struct模块。相反，使用Numpy的结构化数据类型和fromfile。可以查看这里了解更多信息：http://scipy-lectures.github.com/advanced/advanced_numpy/index.html#example-reading-wav-files

• 你可以通过给fromfile传入合适的count=参数，一次性读取所有记录。

像这样（虽然没有测试过，但你能理解这个意思）：

import numpy as np

file = open(input_file, 'rb')
header = file.read(149)

# ... parse the header as you did ...

record_dtype = np.dtype([
    ('timestamp', '<i4'), 
    ('samples', '<i2', (sample_rate, 4))
])

data = np.fromfile(file, dtype=record_dtype, count=number_of_records)
# NB: count can be omitted -- it just reads the whole file then

time_series = data['timestamp']
t_series = data['samples'][:,:,0].ravel()
x_series = data['samples'][:,:,1].ravel()
y_series = data['samples'][:,:,2].ravel()
z_series = data['samples'][:,:,3].ravel()