#!/usr/bin/python
#
# image2wav.py
#
# duncan@linuxbandwagon.com
#
version = '0.1'
# image2wav.py
#
#
# developed on Redhat Linux 7.3 using python v 1.5.2
#
# requires the Python Imageing Library to be installed.
#
# simple program to take an image file and create
# a wav file that will display the image (or an
# approximation) when played through a spectrum analyser.
#
# creates mono 44100 16 bit samples
#
# TODO: selectable output sample length
# more error checking
# map colour images to stereo WAV files
#
# Despite the triviality of this program, it is copyright the author
# under the terms of the Gnu Public License. See
# http://www.gnu.org/licenses/gpl
#
import Image, wave, array, sys, struct
#############################################################################
#
# FFT code from http://www.python.org/topics/scicomp/recipes_in_python.html
#
"""
Find 2^n that is equal to or greater than.
"""
def nextpow2(i):
n = 2
while n < i: n = n * 2
return n
"""
Return bit-reversed list, whose length is assumed to be 2^n:
eg. 0111 <--> 1110 for N=2^4.
"""
def bitrev(x):
N, x = len(x), x[:]
if N != nextpow2(N): raise ValueError, 'N is not power of 2'
for i in range(N):
k, b, a = 0, N>>1, 1
while b >= a:
if b & i: k = k | a
if a & i: k = k | b
b, a = b>>1, a<<1
if i < k: # important not to swap back
x[i], x[k] = x[k], x[i]
return x
"""
FFT using Cooley-Tukey algorithm where N = 2^n. The choice of
e^{-j2\pi/N} or e^{j2\pi/N} is made by 'sign=-1' or 'sign=1'
respectively. Since I prefer Engineering convention, I chose
'sign=-1' as the default.
FFT is performed as follows:
1. bit-reverse the array.
2. partition the data into group of m = 2, 4, 8, ..., N data points.
3. for each group with m data points,
1. divide into upper half (section A) and lower half (section B),
each containing m/2 data points.
2. divide unit circle by m.
3. apply "butterfly" operation
|a| = |1 w||a| or a, b = a+w*b, a-w*b
|b| |1 -w||b|
where a and b are data points of section A and B starting from
the top of each section, and w is data points along the unit
circle starting from z = 1+0j.
FFT ends after applying "butterfly" operation on the entire data array
as whole, when m = N.
"""
def fft(x, sign=-1):
from cmath import pi, exp
N, W = len(x), []
for i in range(N): # exp(-j...) is default
W.append(exp(sign * 2j * pi * i / N))
x = bitrev(x)
m = 2
while m <= N:
for s in range(0, N, m):
for i in range(m/2):
n = i * N / m
a, b = s + i, s + i + m/2
x[a], x[b] = x[a] + W[n % N] * x[b], x[a] - W[n % N] * x[b]
m = m * 2
return x
"""
Inverse FFT with normalization by N, so that x == ifft(fft(x)) within
round-off errors.
"""
def ifft(X):
N, x = len(X), fft(X, sign=1) # e^{j2\pi/N}
for i in range(N):
x[i] = x[i] / float(N)
return x
#######################################################################
#
#
# step 1 - open the image file
#
# some options - mebbe make them command line switchable
# if there is a need.
preserve_aspect = 1
fft_window = 1024
filename = sys.argv[1]
im = Image.open(filename)
print 'Image format and dimensions:'
print im.format, im.size, im.mode
#
# step 2 - get the data as an array that we can feed to
# the fft
# (resize to 1024 wide, for 1024 window thang...)
# convert to grayscale
mim = im.convert('L')
# resize so width is 1/2 of fft window (we pad out the other half with 0s)
if preserve_aspect == 1:
ratio = (fft_window / 2.0) / mim.size[0]
mim = mim.resize((fft_window / 2,int(mim.size[1] * ratio)))
else:
mim = mim.resize((fft_window / 2,mim.size[1]))
pad = [0] * (fft_window / 2)
#
# step 3 IFFT it
#
results = []
print 'transforming.....'
for i in range(mim.size[1]):
print i,' of ',mim.size[1]
slice = []
for j in range(mim.size[0]):
slice.append(mim.getpixel([j,i]))
result = ifft(slice + pad)
results.append(result)
samplelength = len(results[0]) * len(results)
print 'sample length in seconds: ',samplelength / 44100.0
#
# step 4 - save the resultant wav
#
wfile = wave.open(filename+'.wav','w')
wfile.setparams((1,2,44100,samplelength,'NONE',''))
max = 0.0
min = 0.0
print 'finding max amplitude'
for slice in results:
for i in slice:
if abs(i.real) > max:
max = abs(i.real)
print 'max amplitude is ',max
bigint = 2**15
written = 0
print 'writing....'
for slice in results:
print 'written ',written,' of ',samplelength
buff = ''
for i in slice:
val = i.real
buff = buff + struct.pack('h',int(val * bigint / max))
wfile.writeframes(buff)
written = written + (len(buff) / 2)
wfile.close()