wavelet transform on Haskell

久々の更新です。大学卒業したらFFTW扱う暇がなくなって、FFTWモジュールの話は作りかけとまってます。もうしばし待ってください。

で、ずっと更新しないと管理を放棄したかのように思われるので小ネタを。
wavelet変換と呼ばれる信号処理？画像処理？用の変換があります。フーリエ変換と比較して、時間(not周波数)に関しての変換が可能だったりで、多重解像度解析を利用したさまざまな応用があります。

で、普通は信号処理なのでそれっぽい言語を使って実装するところ、Haskellの練習がてら以下のコードを学生時代に書きました。

---------------この区切り表現、Haskellと相性いいかも。
import IO
import System
import List
main = do [filename] <- getArgs
          input <- openFile filename ReadMode
          x <- hGetContents input
          file <- openFile "scaling.dat" WriteMode
-- scaling coefficient
          hPutStr file $ unlines $ map sDouble $
                        [scalingCoefficient level k ((makeList x)!! 0)
                             |k<-[0..((length ((makeList x)!!0)) `div` (2^level))-1 ]]
          hClose file
-- wavelet expansion coefficient
          file <- openFile "wavelet.dat" WriteMode
          hPutStr file $ unlines $ map sDouble $ 
                        [waveletExpansionCoefficient j k ((makeList x)!! 0)
                             | j<-[1..level],
                               k<-[0..((length ((makeList x)!!0)) `div` (2^j))-1 ]]
          hClose file
          hClose input
level :: Int
level = 3
p :: [Double]
p = [0.4829629131445341,0.8365163037378077,0.2241438680420134,-0.1294095225512603]

q :: [Double]
q = zipWith (*) (cycle [1,-1]) (reverse p)
n = 4


scalingCoefficient j k x = if j==0 
                           then (x!!k)
                           else 0.5 * (sum $ zipWith (*) p [scalingCoefficient (j-1) i x| i<-drop (2*k) $cycle [0..((length x)`div`(2^(j-1))-1)]])

waveletExpansionCoefficient j k x = sum $ zipWith (*) q [scalingCoefficient (j-1) i x|
                                                          i<-drop (2*k-1) $cycle [0..((length x)-1)]]

makeList :: String -> [[Double]]
makeList x = transpose $map (map rDouble) $ map words $ lines x

rDouble :: String -> Double
rDouble = read

sDouble :: Double -> String
sDouble = show

---------------

インプットファイルは変換したい信号(波形)を等間隔でサンプリングした結果を、
0.11
0.14
0.15
......
と改行区切りで延々いれたものです。

コンパイルした後
% ghc -o wavele wavelet.hs
% ./wavelet inputfile.name
として使用してください。
scaling係数がscaling.datに、wavelet係数がwavelet.datに出力されます。
ちなみにDaubechies waveletの4。レベル3まで変換します。
信号長さが十分な場合はもっと深いレベルまで変換するのも可。

このコード、作ったのはいいものの信号量が膨大(GByteオーダー)になりすぎてGHCの速度じゃ太刀打ちできなくなって、C++で作り直したんだよなぁ……

それ以前に、汚いコードなのでまねはしないでください。