可以将44.1kHz单通道.wav文件中的一秒读取到长度为44100的数组（称为b）中。给定矩阵A，我们寻求系统Ax = b的解。通过Gauss-Seidel的迭代，向量如果我们将b记录的录音，则将一些白噪声作为我们的初始猜测，并在每次交替中写出Ax，我们会观察到b中高音调的音符首先变得可听，而同时白噪声的音调分解。

最初的12秒.wav文件的音频（白噪声）initialAx.wav

初始Ax，残差和残差FFT的图：

经过一轮迭代，高音转化gauss_seidel_out000000.wav

在光谱中可以看到一些结构：

第二次迭代：gauss_seidel_out000001.wav

第三次迭代：gauss_seidel_out000002.wav

第四次迭代：gauss_seidel_out000003.wav

这一切都在python中完成。将.wav文件加载到数组中，在scipy中还不错。为了避免缓存问题，必须使用稀疏矩阵类，因为12秒的.wav文件需要一个大小为12 * 44100的数组。这是我使用的TridiagonalMatrix类代码片段：

from numpy import *#a tridiagonal matrix class
class TridiagonalMatrix:#initialize with 3 numpy arraysdef __init__(self, upper_in, diag_in, lower_in):self.upper  = upper_inself.diag   = diag_inself.lower  = lower_inself.dim    = diag_in.shape[0]#matrix mulitplicationdef apply(self, v):out = ndarray(self.dim)try:out[0] = self.diag[0]*v[0] + self.upper[0]*v[1]out[self.dim-1] = self.lower[self.dim-2]*v[self.dim-2] + self.diag[self.dim-1]*v[self.dim-1]for i in range(1, self.dim-1):out[i] = self.lower[i-1]*v[i-1] + self.diag[i]*v[i] + self.upper[i]*v[i+1]except(IndexError):print "Wrong sizes"return out

这是处理读取/写入.wav文件然后使用Gauss-Seidel转换为线性系统的代码片段。

from TridiagonalMatrix import *
from numpy import *
from scipy.io import wavfile
import scipy.fftpack
import pylab
import sys
import osdef musical_gauss_seidel(A, b, x0, tol):
"""
do the gauss seidel iteration
but output some sound every now and then..
A is some matrix that lets gauss seidel work
b is a vector that represents a .wav file of a pretty song
x0 is our initial guess for the gauss seidel method (probably random static)
we are going to output the .wav data corresponding to Ax
as Ax gets closer to b (ie the residual gets smaller)
we should hear the song emerge from the initial guess
"""#make noise of the initial approximation to bwavfile.write("gauss_seidel_out000000.wav", 44100, (A.apply(x0)).astype(int16))residual  = A.apply(x0) - b