40年前の卒論の課題をPythonでやってみた

「１／f ノイズの数値解析」　というテーマだったと思います。

　ある時系列データ x[n] を与え、そのパワースペクトルが「１／f 特性」を示すことを数値計算によって示すものです。着手した１９８４年当時、PCは９８０１が使えましたが、言語はBASICのみ。その数年後にはCが普通に使えるようになるのですが、まだBASIC時代に取り組んだ課題です。

　与える時系列データは、"高橋陽一郎先生"（当時東大）が発表されたもので、そのパワースペクトルは１／f 特性を示すことが証明された論文を指導官からもらったのですが学部４年の私にはチンプンカンプン。当時読んだ論文（報告）を入手したいのですが見つかっておりません。すべて手書きで書かれていました。

　その「１／f ゆらぎを生成する、ある関数 f(x) 」は克明に記憶しており、以下のとおりです

　　　　　　f(x) = x  + 2 * x*x ;  (0< x < 0.5),   = 2 *(1 - x ) ;  0.5 ≦ x < 1

 　　　　　　x[1] = f (x[0] ) 　...... 　 x[n+1] = f (x[n] )  のようにして波形を生成します。

このx[n] を時系列波形とみると、そのパワースペクトルが１／f になるというのです。波形を描くとわかるのですが、しばらくゼロに張り付く区間があり時間をかけて漸増し再度動きが始まり・・・を不規則に繰り返します。波形は、「ある初期値」から始まりますが、初期値がわずかでも変化すると波形はまったく別物になります。カオスに見られる性質です。

（まずプログラム制作）

　当時はFFTは簡単にできるものではなく、まず、FFTのアルゴリズム（FORTAN)をBASICで忠実に移植しました。計算時間は、512サンプルで 1分弱だった記憶。しかし試行錯誤をしているとこれでは仕事にならないと思い、ASMで書き換えることにしました。使ったのはPC-9801F2

　当時、PC98のインタプリタ内の「浮動小数点四則演算」のエントリアドレスと使用法が公開されており、その仕様に基づいてASM化。しかし、8086 アセンブラなど手元になく、コーディングシートでプログラムし、ハンドアセンブル、そのコードを　N88-BASIC(86) でDATA文として並べ、poke文でダイレクトにメモリに書き込み、それをBASICからコールする形。たすき掛け演算で使う「関数表」は事前にテーブルを作成しておく方式としました。

　デバッグは・・・MON コマンドなるものが使用できました。原始的なデバッガが装備されています。STEP実行・BREAK・レジスタ表示／変更・メモリダンプ・簡易逆アセンブラ　など。・・・いろいろ苦労したですが、終わりごろはメモリに16進値で直接パッチをあてたり、HEXダンプを見て「脳内逆アセンブル」の域に達したですが。。。

　なんだかんだで、一連のプログラム制作は秋頃までかかり、それから、課題の本計算にはいります。　1000万サンプルを超えるデータを平均化しながら進めるのですが、巨大な配列をつくるにはメモリが足らず、一時的にBASICからデータをFDDに吐き出したこともありました。

　さて計算時間ですが、１つのパラメータセットの計算が大体一昼夜くらい、夜RUNして翌朝に結果がでているくらいの時間感覚です。バグっていたらやり直し orz。これの繰り返し。

　・・・バグもとれて、順調にデータ量産にはいったのが１月。卒論提出が２月中旬で、当然ながら連日徹夜。出力のパワースペクトルグラフは、初期は数値から手書きでグラフを作成していましたが、PC9801のグラフィック＆ハードコピーを作り、最後はひたすら走らす・・・作業にはいっていました。

　記憶は定かでないのですが、CRTに描画さえたグラフを、プリンタにハードコピーする処理も制作したかもしれません。当時のプリンタは「ドット・インパクト方式」で、縦１６ドット・横１２０文字分を一気に印刷でき、専用の制御シーケンスがあった記憶があります。

（今回の検証）

　昭和５９年・６０年ごろは、そんな苦労があったですが、今回、Python 3.9 を使って同じことを再現してみました。

　経過は省略しますが、プログラミングは３時間、実行は２０秒ほどで目的グラフが出力できました。当時使えた主記憶は64kBくらいでしたが、今回は４G実装。CPUクロックに至っては、当時8MHz, 現在1.6GHz という超贅沢。当然現在は、すべてONメモリで計算できます。

（考察）　

　こうやって進化した技術・環境がまた次の成果につながるのだと実感。当時苦労したFFTも、現在ではimport 文にて即利用可能です。私の数年上の先輩は、IBMの大型機で出力した数値が羅列された「ストックフォーム」と格闘しデータをグラフ化していました。

　私の課題では、出たばかりのPC98が使え先輩たちの苦労からみると、まだまだましな方でした。数年後はLatice Cコンパイラが使えるようになり実行速度は飛躍的に恒常。後輩たちはさらにデータを量産しています。私たちはそんな隙間にいた世代でした。

（Pythonプログラム）

今回の試作を公開します。Pythonは得意でないので細部は勘弁してください。

• 実行前に　soundfile を pip にてインストールしておいてください。あとは標準で動作します
• FFTは1024 sampleで処理、10000回実行し「平均スペクトル」を求めます。
• 結果は log-log グラフで表示され、傾きは（ほぼ）−１　となっており、１／f 特性が計算されました。
• 音として聴くため、WAVファイルを生成します。関数の性質上「無音区間」が結構あります。

---

#
# 1-dim 1/f Noise by Takahashi-Map (1984)
#

import soundfile as sf
import numpy as np 
from  matplotlib import pyplot as plt

# Execution
N = 1024        # FFT sample number
Exec = 10000

# Width 
RANGE = range(0, N)
RANGE_s = range(0, int(N/2))

# initital value 
xtmp = 0.10

# Wave array
Wave=[]
Audio =[]

# Spectrum  powerSpectra[n][f_bin]
powerSpectra = []

def saveWave():
    sr = 44100
    filepath = "Takahashi.wav" 
    sf.write(filepath, Audio, sr)

def takahashi(x):
    if (x < 0.5):
        return x + 2 * x * x 
    else:
        return (1-x)*2.0

def GenWave(x0):
    x = x0
    for i in RANGE:
        x = takahashi(x)
        Wave.append(x)  # cleared on each exec.
        Audio.append(x) # add all the samples.
        
    return x
        

# === MAIN =====

for n in range(0,Exec):
    
    print( "- executing {}/{}.".format(n,Exec) )
    
    Wave = []
    xtmp = GenWave( xtmp )

    if False:
        # Graph #
        fig, ax = plt.subplots()

        ax.plot(RANGE, Wave)

        ax.set_xlabel("Time [s]")
        ax.set_ylabel("Signal")
        ax.grid()
        plt.show()

    #
    #  FFT    numpy.fft.fft(a, n=None, axis=-1, norm=None)
    #

    Sp = np.fft.fft( Wave, n=N, axis=-1, norm=None)

    power = (np.abs(Sp))**2 / N 
    powerSpectra.append(power)

    if False:
        plt.plot(power)
        plt.show()

#
# Averaging
#

print( "- averaging.")
aveSpectra = []
for bin in RANGE_s:
    p = 0
    for n in range(0,Exec):
        p = p + powerSpectra[n][bin]
    
    p = p / Exec
    aveSpectra.append(p)
#

# create WAV
saveWave()

# averaged spectrum
if True:
    
    plt.figure(figsize=(6,4))
    plt.title("1/f spectrum by Takahashi-map")
    plt.xscale('log')
    plt.yscale('log')
    
    plt.grid(which='major',color='black',linestyle='-')
    plt.grid(which='minor',color='gray',linestyle='--')
 
    plt.plot(aveSpectra[0: int(N/2)])
    plt.show()

---