From Karaoke

女房に誘われて、久しぶりに からおけ へ行ってきた。入場料がわりに頼む必要 がある、ドリンク(430円)と、一人80円/時間の部屋代だけで済むとは、一体これで 儲かるのかと心配になる。

一曲毎に採点してしてくれるってのは、昔有ったかなあ？ 今日の最高得点は、89 点だった。しつこく同じ曲を何度も歌ったけど、採点はぶれなかった。一体、どう いう仕組みなんだろう？ ちょっと気になったぞ。後、曲毎に消費カロリーが 表示されてたけど、こちらは、ぶれが無い事からすると、DBにでも保存してるのを ただ表示してるだけだな。カロリーを気にするおばさん対策と思うぞ。

ちなみに、最高得点だった曲は、「鉄腕アトム」。この曲には思い出があって、 昔、東南アジアからの研修生受け入れの宴会で歌った時、大受けして彼らと 何度も歌ったんです。中国とか台湾にもしっかりアニメが浸透してるなあ。

と言う事で、少々音関係に興味が出てきた。マイクで音を拾って、波形を見る。 その波形を、周波数ドメインで見るなんて事をやってみたい。

波形の作成

が、マイクは、生憎持ち合わせていないので、取り合えずサイン波を合成してやればいいだろう。 どうせやるなら、昔学校で習った、AM放送の原理を再現してみる。 さて、言語は何にしよう？ 最近触っていない ruby でやってみる。

#!/usr/local/bin/ruby
DATAN = 256         # Wave size 

Sin_amp = {  # Table of Freq vs Amp,   freq => amp
     3 => 0.3,  
}

Car = 64

def cal_point(n)
  k =0.0
  Sin_amp.each do |freq, amp| f = freq.to_i
     k = 1 + amp * Math::sin(2 * Math::PI * f * n / DATAN)
  end
  val = k * Math::cos(2 * Math::PI * Car * n / DATAN)
  val
end

DATAN.times do |n|  
  printf("%.6f\n", cal_point(n) )
end

TBSなら954KHzだけど、このプログラムでは64（と言う、搬送波）、アナウンサーの声 の変わりに、周波数３と言う信号波、振幅は、搬送波に対して0.3倍（だから、変調率で言うと、 0.3で、良かったかな？ 中学生の時、アマチュア無線の資格を取った時以来だから 、ちょっと自信なさげ）

出力は、電波じゃなくて、標準出力へ放出される。

FFT するぞ

次は、スペクトラムアナライザーの心臓部、FFTだ。うーん。FFTねぇ。ネットを 探してみると、rubyだけで書かれた例も見つかったけど、上と同じRubyを使うのも なんだかつまんないので、別の路線を考える。

思い出したのが、「プログラミング言語SCHEME」に出ていた例 これを頂いてきて、読み込みルーチンを付け加えてやろう。最近ご無沙汰してる gaucheを使う事にする。

#!/usr/local/bin/gosh
;;; Original by fft.ss
;;; Copyright (C) 1996 R. Kent Dybvig
;;; from "The Scheme Programming Language, 2ed" by R. Kent Dybvig

(define DATA_FILE "wavein.txt")

(define (w-powers n)
  (let ((pi (* (acos 0.0) 2)))
    (let ((delta (/ (* -2.0i pi) n)))
      (let f ((n n) (x 0.0))
	(if (= n 0)
	    '()
	    (cons (exp x) (f (- n 2) (+ x delta))))))))

(define (evens w)
  (if (null? w)
      '()
      (cons (car w) (evens (cddr w)))))

(define (interlace x y)
  (if (null? x)
      '()
      (cons (car x) (cons (car y) (interlace (cdr x) (cdr y))))))

(define (split ls)
  (let split ((fast ls) (slow ls))
    (if (null? fast)
	(values '() slow)
	(call-with-values
            (lambda () (split (cddr fast) (cdr slow)))
	  (lambda (front back)
	    (values (cons (car slow) front) back))))))

(define (butterfly x w)
  (call-with-values
      (lambda () (split x))
    (lambda (front back)
      (values
       (map + front back)
       (map * (map - front back) w)))))

(define (rfft x w)
  (if (null? (cddr x))
      (let ((x0 (car x)) (x1 (cadr x)))
        (list (+ x0 x1) (- x0 x1)))
      (call-with-values
	  (lambda () (butterfly x w))
        (lambda (front back)
          (let ((w (evens w)))
            (interlace (rfft front w) (rfft back w)))))))

(define (fft x)
  (rfft x (w-powers (length x))))


;;; below line's add by S.K. ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

(define (make-list file)
  (define (build lis x fp)
    (if (eof-object? x)
	(reverse lis)
	(build (cons x lis) (read fp) fp)))
  (call-with-input-file file
    (lambda (fp)
      (build '() (read fp) fp))))

(define (other-way-make-list file)  ;; Care order of val list
  (with-input-from-file file
    (lambda ()
       (port-fold-right (lambda (x lis) (cons x lis)) '() read))))

(define (do-fft file)
  (let ([x (make-list file)])
    (for-each (lambda (v) (format #t "~s\n" (abs v)))
	       (fft x))))

(do-fft DATA_FILE)

入力は、ファイル渡しで、出力は、標準出力に出す事にした。

表示部の作成

次は、結果の表示だ。グラフを書く事になるので、スマートに済ませてしまおう。 Rを使って、一機にpngファイルを作る事にする。

#!/usr/local/bin/Rscript

# pdf(file="Rplots.pdf")
png(file="Rplots.png")

par(bg="grey95", mar=c(4,4,1,1), mfrow=c(2,1))

wave <- read.table("wavein.txt")
plot(wave$V1, type="l")

spec <- read.table("specout.txt")
n <- length(spec$V1)
s  <- spec$V1[1:(n/2)]
plot(s/n, type="h")

dev.off()

お好みで、PDFにも吐き出せるぞ。

全部をまとめてみる。

さて、部品は揃った。後は、これらを結合するんだけど、相場は shell を使う んだろうけど、後々の事（拡張）を考えて、Makefileにしておく。 そうすれば、make am とか（今回）、FM波形が必要なら、適当に波形を合成して make fm とするだけだ。

am:	
	./AMgen.rb    > wavein.txt
	./FFT.scm     > specout.txt
	./disp.R
	xv Rplots.png
	#acroread Rplots.pdf

結果発表

某メーカーもびっくりの、よくばり表示器です。

オシロFFT

参考にデジタル信号処理 を上げておきます。