#!/usr/local/bin/wish -f
#
#  Fourier Transform Training Software
#
#  coded by Shigeki Sagayama, JAIST, 20 Sep 1998
#
#  今後の追加内容案: log spectrum, 信号値のエディット (text widget)？
#

# set Date "20 Sep 1998"
set Date "24 Oct 1998"
puts stdout \
"< Fourier Transform Training Software (S. Sagayama, JAIST, $Date) >"
#------------------------------------------------------------------------------
# help
#------------------------------------------------------------------------------
proc helpIntro {} { helpDisp { Help: Introduction - FourierTrans の概要

                北陸先端科学技術大学院大学  嵯峨山 茂樹

これは、北陸先端科学技術大学院大学の基幹科目講義 I212 情報解析学特論の受講者
を対象に、フーリエ変換を復習するための自習用ソフトウェアとして作成したもので
す。

フーリエ変換の性質を感覚的に把握するためのもので、与えられた信号波形に対する
フーリエスペクトルを表示します。ただし、これは DFT (discrete Fourier
transform, 離散フーリエ変換) の計算によっているので、細かい点では定義通りの
フーリエ積分変換とは若干の差があることは了解しておいて下さい。

このプログラムは Tcl/Tk のスクリプトとして書いてあります。機種やＯＳに大きく依
存せず、普及している無償のソフトのみを用いて、画面表示に動きのある自習プログラ
ムを作るには Tcl/Tk が適当であると考えたからです。しかし Tcl/Tk は数値計算や信
号処理などへの応用を想定していないので、実行速度が遅いのは仕方ありません。

作者は、教材作りを目的に Tcl/Tk を勉強し始めてから二週間ほどの時にこのプログ
ラムを書いたので (後に若干の修正を行なっていますが)、プログラムの書き方が未
熟な箇所は随所に見られることと思います。しかし、作者に無断で改善・改造して、
再配布することは避けて下さい。再配布に関して詳しくは Help: Copyright をクリッ
クして見て下さい。

以上。}
}
#------------------------------------------------------------------------------
proc helpUsage {} { helpDisp {\
Help: Usage - FourierTrans の使用法

使い方は簡単です。ほとんどマウス操作だけで済みます。

(1) 画面の説明
==============
    およそ下半分の画面は時間領域、すなわち解析すべき信号の画面です。お
    よそ上半分の画面は、周波数領域、すなわち信号のフーリエ変換です。
    時間領域を表の世界、フーリエ変換を裏の世界にたとえます。また、ポジ
    とネガの関係にたとえます。そのため、このソフトでも周波数領域を黒地
    にしました。

    時間軸も周波数軸も、中央が0に相当します。つまり、左半分は負の時刻
    あるいは周波数、右半分は正の時刻あるいは周波数です。

    下の方に、スペクトル表示のスケーリングと時間のスケーリングを変化さ
    せるためのスケール窓があります。

    最下段には各種のボタンスイッチ類が配置されています。

(2) 信号入力
============
    任意の信号波形のフーリエ変換を観察できます。信号波形の入力の仕方は
    次のように4通りあります。

    1: マウスによりサンプル値を設定

      時間領域画面の緑の小さい正方形(128個ある)を、マウスボタン１を押
      しながらひきずって(drag する)上下方向に動かして行ないます。入力
      中のサンプルは色が赤くなり、その値は中央の x(t)=... に表示されま
      す。

    2: ファイルから読み込み

      時間軸の順に一行ごとにサンプル値が ASCII コードで書かれた128行の
      ファイルから入力するものです。"load data" ボタンをクリックすると
      ファイル名を尋ねてくるので、まずファイル名窓でマウスボタン1をク
      リックしてからキーボードからファイル名を投入し、リターンキーを押
      します。一行ずつ処理が進んで行くのが見られます。ファイル名投入を
      キャンセルするにはマウスボタン3をクリックします。

    3: 信号関数を定義 ("def func" をクリック)

      ユーザ関数の式を与えて、それを信号として入力するものです。"def
      func" ボタンをクリックすると関数式を尋ねてくるので、まずマウスボ
      タン1をクリックしてから、キーボードから関数式を投入し、リターン
      キーを押します。一行ずつ処理が進んで行くのが見られます。ファイル
      名投入をキャンセルするにはマウスボタン3をクリックします。数式の
      与え方は help の別項目を見て下さい。要点は時刻を $t として与える
      ことです。(例: 100*sin(3.14*$t) )

    4: クリアする ("clear" をクリック)

      一挙に全零入力(最初の状態)に戻したい場合は、"clear" ボタンをクリッ
      クします。

(3) 信号表示
============
    上のような方法で入力できる信号波形の表示について、以下のようなこと
    が可能です。

    1: 時間スケールの調整 ("time scale" スケールを動かす)

      信号波形に縦に引かれている時間軸の目盛間隔を変化させることができ
      ます。周波数軸の目盛も自動的に変化します。たとえば時間目盛を mS 
      単位と考えれば、周波数目盛の単位は kHz となる。時間スケールを調
      整することで、信号波形形状はそのままで時間長さを変えてそのフーリ
      エ変換を考えます。

    2: 入力用の小さい緑の正方形を消す ("hide pts" をクリック)

      関数形をよく観察したり、印刷するなどの目的で、入力用の小さい緑
      の正方形が邪魔な場合は、"hide pts" をクリックをします。もう一度
      同じボタン("show pts" となっているはず)をクリックすると、表示が
      復活します。もちろん、表示消去中は、マウスからの信号入力はでき
      ませんが、他の3種の入力方法は有効です。

    3: 滑らかな曲線表示 ("smooth" ボタンをクリック)

      "smooth" ボタンをクリックすると、信号波形もスペクトル形状も滑ら
      かな曲線 (Bezier 曲線) で表示されます。表示の滑らかさが必要な時
      に有用ですが、サンプル値が正確に表現されないことに注意する必要
      があります。もう一度同じボタン("polyln" となっているはず)をクリッ
      クすると、もとの折れ線による表示モードに戻ります。

    4: eps ファイル作成 ("save eps" ボタンをクリック)

      信号の形の画面の eps (encapsulated PostScript) ファイルを作成す
      るものです。印刷に、あるいは TeX 文書への埋め込みなどの用途が考
      えられます。下段の "save eps" ボタンをクリックすると eps ファイ
      ル名を尋ねてくるので、まずマウスボタン1をクリックしてから、キー
      ボードからファイル名を投入し、リターンキーを押します。ファイル名
      投入をキャンセルするにはマウスボタン3をクリックします。

(3) フーリエ変換表示
====================
    周波数領域(フーリエ変換)は、信号の入力に伴って逐次表示が変化します。
    赤い曲線は実数部、青い曲線は虚数部を表します。この表示に関して以下
    のことが可能です。

    1: スペクトル値スケールの調整 ("spectrum scaling factor" スケール
       を動かす)

      "spectrum scaling factor" スケールを動かすと、スペクトル値をスケー
      リング値で割って表示します。見やすい値に調整して下さい。

    2: 周波数スケールの調整 ("time scale" スケールを動かす)

      "time scale" スケールを動かすと、時間軸の目盛間隔と同時に周波数
      の目盛間隔も変化します。

    3: 絶対値の表示 (黄色)

       "show abs" ボタンをクリックすると、フーリエ変換の絶対値(実数部
       と虚数部の二乗和の平方根)、すなわちパワースペクトルの平方根が、
       黄色で表示されます。同じボタンをもう一度押すと絶対値表示は消え
       ます。絶対値表示をしないほうが、動作が若干速くなります。

    4: 滑らかな曲線表示

       "smooth" ボタンをクリックすると、信号波形もスペクトル形状も滑ら
       かな曲線 (Bezier 曲線) で表示されます。表示の滑らかさが必要な時
       に有用ですが、サンプル値が正確に表現されないことに注意する必要
       がありますが、スペクトルの表示の場合は有用なことが多いようです。

    5: eps ファイル作成 ("save eps" ボタンをクリック)

      スペクトル画面の eps (encapsulated PostScript) ファイルを作成す
      るものです。印刷に、あるいは TeX 文書への埋め込みなどの用途が考
      えられます。"save eps" ボタンをクリックすると eps ファイル名を尋
      ねてくるので、まずマウスボタン1をクリックしてから、キーボードか
      らファイル名を投入し、リターンキーを押します。ファイル名投入をキャ
      ンセルするにはマウスボタン3をクリックします。

以上。}
}
#------------------------------------------------------------------------------
proc helpCopyright {} { helpDisp {\
Help: Copyright

◆ このプログラムの著作権は、作者・嵯峨山茂樹 (sagayama@jaist.ac.jp) に属し
ます。

このプログラムは、北陸先端科学技術大学院大学の専門科目講義 I212 情報解析学特
論の受講者を対象に、フーリエ変換を理解するための自習用ソフトウェアとして嵯峨
山が作成したものです。

作者は、教材作りを目的に Tcl/Tk を勉強し始めて、二週間ほどでこれを含む数個の
プログラムを書いたので、プログラムの書き方が下手な箇所は随所に見られることと
思います。また、フーリエ変換のような信号処理分野を Tcl/Tk のスクリプトで行な
うという発想がそもそも特異なので、改善の余地は多くあることと思います。皆様の
暖かい改善の御意見を歓迎します。

このような不完全な代物ではありますが、しかし、作者に無断で改善・改造して、他
の人に配布することは避けて下さい。そこで、以下に使用条件を述べます。

◆ このソフトウェアを改変することなくそのまま全部を再配布することは許可しま
す。大いに使って頂いて結構です。ただし、改善・改変して再配布しないで下さい。
この著作権に関する宣言の部分も含め、一字一句変えずにこのままならば再配布は自
由です。但し、漢字コードの変換などは構いません。

◆ 作者の許可なく改変(一部分のみの使用を含む)したものは、再配布することを禁
じます。

◆ バグやプログラミングに関する御意見は、作者に御連絡下さい。バグ情報、改善
意見、改良の努力、利用経験、応用情報などは歓迎致します。

以上。}
}
#------------------------------------------------------------------------------
proc helpUserfunc {} { helpDisp {\
help: ユーザ定義関数

◆ "def func" と表示されたボタンをクリックすると、ユーザ定義関数を尋ねてきま
す。時刻は $t で表されます。そこで、Tcl の expr コマンドの文法に従って、式を
与えて下さい。

◆ 例えば、

    cos 関数:      100*cos(10*$t)

です。係数の 100 は、画面の中での表示の都合です。また、

    sinc 関数:      100*(sin(10*$t)+0.01)/(10*$t+0.01)

とするとよいでしょう。式中の 0.01 は、分母・分子がともに 0 になるのを防ぐた
めです。もちろん場合分けでも可能です。

◆ Tcl の数式は、C の文法に近く、かなり自由に表現できます。一部を紹介すると、

   算術演算    -a, a*b, a/b, a%b, a+b, a-b

   数学関数    abs(x), acos(x), asin(x), atan(x), atan2(x,y), ceil(x),
               cos(x), soch(x), double(i), exp(x), floor(x), fmod(x,y),
               hypot(x,y), int(x), log(x), log10(x), pow(x,y), round(x),
               sin(x), sinh(x), sqrt(x), tan(x), tanh(x)

   論理式      !a,  a<b, a>b, a<=b, a>=b, a==b, a!=b, a&&b, a||b,

   条件分岐    a?b:c (a は論理式、b と c は場合ごとの数式)

   ビット演算  ~a, a<<b, a>>b, a&b, a~b, a|b,

以上。}
}
#------------------------------------------------------------------------------
proc helpSavesig {} { helpDisp {\
help: ファイルへの信号データ出力

    信号画面上の "save data" ボタンをクリックすると、保存出力用のファイル名
を尋ねて来ます。ファイル名投入の窓の中でマウスボタン 1 をクリックすると、キー
ボードからの入力を受け付けるようになります。ファイル名を入れてリターンキーを
押すと、指定されたファイルに、現在表示されている信号のサンプル値系列が保存さ
れます。保存の形式は、ASCII 数字のみで、一行ごとに1データが書き込まれます。

   ファイル名を尋ねてきたところでマウスボタン3をクリックすれば、キャンセルで
きます。

以上。}
}
#------------------------------------------------------------------------------
proc helpEpssig {} { helpDisp {\
help: 信号画面の eps データのファイルへの出力

    信号画面上の "save eps" ボタンをクリックすると、保存出力用のファイル名を
尋ねて来ます。ファイル名投入の窓の中でマウスボタン 1 をクリックすると、キー
ボードからの入力を受け付けるようになります。ファイル名を入れてリターンキーを
押すと、指定されたファイルに、eps ファイルがが保存されます。

   ファイル名を尋ねてきたところでマウスボタン3をクリックすれば、キャンセルで
きます。

以上。}
}
#------------------------------------------------------------------------------
proc helpSavespec {} { helpDisp {\
help: スペクトル画面のデータ出力

    スペクトル画面上の "save data" ボタンをクリックすると、保存出力用のファ
イル名を尋ねて来ます。ファイル名投入の窓の中でマウスボタン 1 をクリックする
と、キーボードからの入力を受け付けるようになります。ファイル名を入れてリター
ンキーを押すと、指定されたファイルに、現在表示されている信号のサンプル値系列
が保存されます。保存の形式は、ASCII 数字のみで、一行ごとに1データが書き込ま
れます。

   ファイル名を尋ねてきたところでマウスボタン3をクリックすれば、キャン
セルできます。

以上。}
}
#------------------------------------------------------------------------------
proc helpEpsspec {} { helpDisp {\
help: スペクトル画面の eps ファイルへの出力

    スペクトル画面上の "save eps" ボタンをクリックすると、保存出力用のファイ
ル名を尋ねて来ます。ファイル名投入の窓の中でマウスボタン 1 をクリックすると、
キーボードからの入力を受け付けるようになります。ファイル名を入れてリターンキー
を押すと、指定されたファイルに、画面の eps ファイルが保存されます。

   ファイル名を尋ねてきたところでマウスボタン3をクリックすれば、キャン
セルできます。

以上。}
}
#------------------------------------------------------------------------------
proc helpLoadsig {} { helpDisp {\
help: ファイルからの信号データ入力

    信号画面上の "load data" ボタンをクリックすると、保存出力用のファイル名
を尋ねて来ます。ファイル名投入の窓の中でマウスボタン 1 をクリックすると、キー
ボードからの入力を受け付けるようになります。ファイル名を入れてリターンキーを
押すと、指定されたファイルから信号のサンプル値系列を読み込み、フーリエ変換を
計算して表示します。ファイルの形式は、ASCII 数字のみで一行ごとに1データが書
き込まれている形式です。"save data" で得られた退避データ形式がそのまま使われ
ます。

   ファイル名を尋ねてきたところでマウスボタン3をクリックすれば、キャンセルで
きます。

以上。}
}
#------------------------------------------------------------------------------
proc helpSaveLoad {} { helpDisp {\
help: ファイル入出力

ファイルとのやりとりは、ごく簡単な以下の種類があります。help の操作説明の部
分にも説明してありますが、まとめると、

(1) ファイルへの信号データ出力

    信号画面上の "save data" ボタンをクリックすると、保存出力用のファイル名
を尋ねて来ます。ファイル名投入の窓の中でマウスボタン 1 をクリックすると、キー
ボードからの入力を受け付けるようになります。ファイル名を入れてリターンキーを
押すと、指定されたファイルに、現在表示されている信号のサンプル値系列が保存さ
れます。保存の形式は、ASCII 数字のみで、一行ごとに1データが書き込まれます。

   ファイル名を尋ねてきたところでマウスボタン3をクリックすれば、キャンセルで
きます。

(2) ファイルからの信号データ入力

    信号画面上の "load data" ボタンをクリックすると、保存出力用のファイル名
を尋ねて来ます。ファイル名投入の窓の中でマウスボタン 1 をクリックすると、キー
ボードからの入力を受け付けるようになります。ファイル名を入れてリターンキーを
押すと、指定されたファイルから信号のサンプル値系列を読み込み、フーリエ変換を
計算して表示します。ファイルの形式は、ASCII 数字のみで一行ごとに1データが書
き込まれている形式です。"save data" で得られた退避データ形式がそのまま使われ
ます。

   ファイル名を尋ねてきたところでマウスボタン3をクリックすれば、キャンセルで
きます。

(3) 信号画面の eps データのファイルへの出力

    信号画面上の "save eps" ボタンをクリックすると、保存出力用のファイル名を
尋ねて来ます。ファイル名投入の窓の中でマウスボタン 1 をクリックすると、キー
ボードからの入力を受け付けるようになります。ファイル名を入れてリターンキーを
押すと、指定されたファイルに、eps ファイルがが保存されます。

   ファイル名を尋ねてきたところでマウスボタン3をクリックすれば、キャンセルで
きます。

(4) スペクトル画面のデータ出力

    スペクトル画面上の "save data" ボタンをクリックすると、保存出力用のファ
イル名を尋ねて来ます。ファイル名投入の窓の中でマウスボタン 1 をクリックする
と、キーボードからの入力を受け付けるようになります。ファイル名を入れてリター
ンキーを押すと、指定されたファイルに、現在表示されている信号のサンプル値系列
が保存されます。保存の形式は、ASCII 数字のみで、一行ごとに1データが書き込ま
れます。

   ファイル名を尋ねてきたところでマウスボタン3をクリックすれば、キャン
セルできます。

(5) スペクトル画面の eps ファイルへの出力

    スペクトル画面上の "save eps" ボタンをクリックすると、保存出力用のファイ
ル名を尋ねて来ます。ファイル名投入の窓の中でマウスボタン 1 をクリックすると、
キーボードからの入力を受け付けるようになります。ファイル名を入れてリターンキー
を押すと、指定されたファイルに、画面の eps ファイルが保存されます。

   ファイル名を尋ねてきたところでマウスボタン3をクリックすれば、キャン
セルできます。

以上。}
}
#------------------------------------------------------------------------------
proc helpLearn {} { helpDisp {\
help: How to Learn - これを用いた学習法のヒント

このソフトウェアを用いてフーリエ変換の学習をしましょう。たとえば、つぎのよう
なフーリエ変換の性質を実際に試してみるとよいでしょう。

(1) 対称関数         フーリエ変換の実数部(桃色)は対称で、虚数部(水色)は全部
                     ０になる。

(2) 反対称関数       フーリエ変換の実数部(桃色)は全部 0 で、虚数部(水色)は反
                     対称になる。

(3) 矩形関数         フーリエ変換は sinc 関数になる。パワースペクトルは sinc
                     関数の２乗。

(4) sinc 関数        フーリエ変換は矩形関数になる。

(5) 等間隔パルス列   フーリエ変換も等間隔パルス列になる。

(6) 時間遅れ         フーリエ変換は変わるが、パワースペクトルは同じ。

(7) 三角パルスのフーリエ変換は？

(8) ガウス波形のフーリエ変換は？

(9) 同じ長さの波形パルスで、スペクトルの広がりが最も小さくできる形は何？

(10) 一般的に、信号の長さとスペクトルの広がりの関係は？

以上。}
}
#------------------------------------------------------------------------------
proc helpFFT {} { helpDisp {\
help: Use FFT - FFT を用いた高速フーリエ変換

Tcl/Tk は計算処理や信号処理には向かない言語です。(フーリエ変換のような信号処
理の自習プログラムを Tcl/Tk で作るということ自体、とんでもない考えであって、
実際、作者はかなり苦労しました。) そのため、本プログラムでのフーリエ変換の計
算は時間がかかります。

入力信号データ、あるいは信号関数定義の場合に、高速フーリエ変換(FFT)の外部プ
ログラムを用いて高速に結果を得るモードを用意してあります。それを実行するため
には、実行する前にプログラム "FFT128" をインストールする必要があります。イン
ストールの仕方は、一度だけメニュー項目 "install" をクリックして数秒待つだけ
です。作業ディレクトリを使って C のソースプログラムをコンパイルし、カレント
ディレクトリに "FFT128" を作成します。以後は、必要に応じて "FFT128" の置き場
所を、環境変数 PATH に含まれる任意のディレクトリに変えても構いません。

以上。}
}
#------------------------------------------------------------------------------
proc helpEnglish {} { helpDisp {\
Help: General information and usage in English

                    Shigeki Sagayama, Professor, JAIST

This is a visualization software for students taking Course I212 ("Information
Analysis") at JAIST.  Change the signal by dragging the small yellow boxes
with mouse button 1.  The resulted Fourier transform will be displayed
immediately.  Note that the Fourier transform displayed here is only an
approximation to the exact integral Fourier transform through discrete Fourier
transform calculations.

The purpose of this software is to provide a tool for intuitive understanding
of Fourier transform.  Be noted that this software adopts DFT (discrete
Fourier transform) and results in a small differences from exact computation
of Fourier integrals.

It is very easy to use.

(1) Screens
===========
    The lower half screen is the time domain, i.e., the signal to analyze.
    The upper half screen is the frequency domain, i.e., the spectrum of the
    given signal which is the Fourier transform of the signal.  As the time
    domain is often compared to photo positives, while the frequency domain
    is to negatives, the screens are positive for the time domain and
    negative for the frequency domain.

    Both time and frequency axes has their origins at the centers of
    screens.  Namely, the left half represents negative time or frequency,
    while the right half represents positive time or frequency.

    Time and frequency scaling can be changed by manupilating the "time
    scale" at the bottom of the screen.

    Buttons provided at the bottom are used for changing display modes, file
    I/O, and defining the user function.

(2) Signal input
================
    The Fourier transform can be observed for an arbitrary input signal
    waveform.  There are four ways of specifiying the signal waveform.

    1: Specify the signal sample sequence using the mouse

      Move as appropriate the small yellow squares (128 squares are shown)
      by dragging with mouse button 1.  The sample point turns red and the
      sample value is displayd as "x(t)=..." in the box at the center.
      
    2: Read a data file (click "load data")

      128 ASCII code values of samples are read sequentially line by line
      from a data file.  Click "load sig" button and the program asks for a
      file name entry.  Activate the entry window by clicking mouse button
      1.  Enter the file name and terminate it by the <return> key.
      Line-by-line input will be seen as Fourier transform calculation going
      on.  In the case of cancelation of file name entry, click the "dismiss"
      button instead.

    3: Define the user function (click "def func")

      Provide a mathematical formula to define the user signal function of
      time $t.  Click the "user func" button and the program will ask for a
      function.  Activate the entry window by clicking mouse button 1.
      Enter the function and terminate it by the <return> key.  Signal
      samples will be generated one by one along with the Fourier transform
      calculation.  In the case of cancelation of the function entry, click
      "dismiss" button instead.  Note that the time is represented by "$t"
      witch must be included in the user function definition.  For example,
                100*sin(3.14*$t)
      Refer to the Tcl math syntax for more information.  (Japanese version
      for the outline is available by choosing help : user function.
      
    4: Clear the data (click "clear")

      Clear the signal by changing all samples into zeros by clicking the
      "clear" button.  The spectrum will be also cleared, consequently,

(3) Signal display
==================
    Signal display has several options as follows:

    1: Adjusting the time scale (move "time scale")

      The time scale can be changed by dragging the nob of "time scale".
      The consequent frequency scale also changes automatically.  For
      example, if the time is in msec, the frequency scale is in kHz.
      Preserving the waveform, changing the time scale gives the frequency
      scaled Fourier transform.

    2: Hide the small squares (click "hide pts")

      Click "hide pts" to remove the small squares overlapping on the signal
      waveform if only the waveform is preffered for observing the shape or
      printing the waveform.  Clicking the same button (appearing with the
      sign of "show pts") will restore the small squares.  While the small
      squares are unseen, the manual input of signal is inhibited, but other
      three modes are still available.

    3: Smooth curve display (click "smooth")

      Clicking the "smooth" button changes smooth the displayed waveform and
      the spectral shapes.  Note that these Bezier curve may not be accurate
      in representing the sample sequences.  Clicking the same button
      (appearing as "polyln") returns to the original polyline mode.

    4: Creating an eps file (click "save eps")

      An eps (encapsulated PostScript) file can be generated from the signal
      waveform.  Applications included printing the waveform and embeddeed
      figures in TeX documents.  Click the "sig eps" button and the program
      will ask for a file name to save the eps code.  Activate the file name
      entry window by clicking mouse button 1 in the window.  Enter the file
      name and terminate it with <return>.  The "dismiss" button cancels the
      file name entry on the way.

(3) Fourier transform display
=============================
    Fourier transform display changes simultaneously along with the
    signal input.  The red line is the real part of the transform and
    the blue line represents the imaginary part.  The following
    options are available for the Fourier transform display.

    1: Scaling the spectrum (move "spectrum scaling factor")

      Moving the scale indicated "spectrum scaling factor" changes the
      spectrum scaling.  Adjust the scale so as to see the spectrum shape
      appropriately.

    2: Scaling the frequency (move "time scale")

      Moving the "time scale" changes the time scale and the frequency scale
      simaultaneously.

    3: Power spectrum display
 
      The square root values of the power spectrum is displayed in a yellow
      line by clicking the "show abs" button.  Clicking the same botton
      (showing a sign "hide abs") removes this absolute value plot.
      Displaying the power spectrum increases the computational load.

    4: Smooth curves

      Clicking the "smooth" button forces the smooth lines (Bezier curves)
      for signal and spctrum lines.  This is useful when smoothness is
      necessary, though the sample values may not be accurate.  Recommended
      for spectrum printing.

    5: Creating an eps file (click "save eps")

      An eps (encapsulated PostScript) file for the displayed spectrum
      shapes can be obtained.  Applications include printing and embbeding a
      figure in TeX documents.  Click the "spec eps" and the program will
      ask for a fiel name to save the eps code.  Activate the window by
      clicking the button 1 and enter a file name terminated by <return>.
      To cancel the entry before typing a <return>, click the "dismiss" botton.

End.}
}
#------------------------------------------------------------------------------
### Create a new message window to show the specified contents ###
proc helpDisp { contents } {
  set w .basic
  catch {destroy $w}
  toplevel $w
#  dpos $w
  wm title $w "Fourier Transform"
  wm iconname $w "Fourier Transform"
  button $w.ok -text OK -command "destroy $w" -fg black -bg yellow
  text $w.t -relief raised -bd 2 -yscrollcommand "$w.s set" -setgrid true \
    -fg black -bg white
  scrollbar $w.s -relief flat -command "$w.t yview"
  pack $w.ok -side bottom -fill x
  pack $w.s -side right -fill y
  pack $w.t -expand yes -fill both
  $w.t insert 0.0 $contents
  $w.t mark set insert 0.0
  bind $w <Any-Enter> "focus $w.t"
}
#------------------------------------------------------------------------------
# initial settings
#------------------------------------------------------------------------------
# constants - can be changed by command arguments
set W 772
set H 320
set N 128

### check command arguments ###
for { set i 0 } { $i < $argc } { incr i } {
  switch [lindex $argv $i] {
    "+length" { incr i; set N [lindex $argv $i] }
    "+width" { incr i; set W [lindex $argv $i] }
    "+height" { incr i; set H [lindex $argv $i] }
    default {
      puts stdout "options: +length L, +width W, +height H"
      exit
    }
} }

set DivT0 60
set Resol [expr floor($W/$N)]
set X0 [expr ($W-($N-1)*$Resol)/2+2]
set Y0 [expr $H/2]
set Xorig [expr $X0+floor($W/$Resol/2)*$Resol]
set Yorig $Y0

# Frequency range [-pi, pi]
set OmegaMax [expr 3.1459265359*2]
#  empirically, 4600.0 was chosen to fit the relation: time x frequency = 1
set TimeFreq 4600.0

# Magnify the central part: frequency range changed into [-pi/2, pi/2]
set OmegaMax [expr $OmegaMax/2]
set TimeFreq [expr $TimeFreq*2]

# computation variables
# initial values - Re, Im, Power Spec
for { set i 0 } { $i < $N } { incr i } {
  set Signal($i) 0
  set SpecRe($i) 0
  set SpecIm($i) 0
  set SpecPw($i) 0
}

set Buff { }
set PrevN 0
set PointList { }

### internal buffer variables
# xy-coord list
set Scale 10
set XYSignal { }
set XYSpecRe { }
set XYSpecIm { }
set ThisY "none"

set Points(0) 0

set CurvSig 0
set CurvRe 0
set CurvIm 0
set ShapeSig { }
set This "none"
set Seen 0
set Smooth 0
set ShowAbs 0
set ShowPnts 1
set AutoScale 1
set GridShown 1
set DivT $DivT0
set UseFFT 0
#------------------------------------------------------------------------------
# display setup
#------------------------------------------------------------------------------
### window title ###
wm title . "Fourier Transform (Sagayama)"

#### canvas ####
canvas .c -width [expr $W+10] -height $H -bg #fff0bc
canvas .d -width [expr $W+10] -height $H -bg #00000f

#------------------------------------------------------------------------------
### buttons: start, stop, quit, time, numb ###
frame .f

label .f.m -relief raised \
  -fg #ffffff -bg #193005 -text \
"Learn Fourier Transform (sagayama@jaist.ac.jp, $Date)"

menubutton .f.help -text "help" -menu .f.help.m -relief raised \
  -fg #ffffff -bg #ff1c35
menu .f.help.m -fg #000000 -bg #ff1c35
.f.help.m add command -label "English" -command helpEnglish
.f.help.m add command -label "Introduction" -command helpIntro
.f.help.m add command -label "How to use it" -command helpUsage
.f.help.m add command -label "Copyright" -command helpCopyright
.f.help.m add command -label "User function" -command helpUserfunc
.f.help.m add command -label "Save/Load" -command helpSaveLoad
.f.help.m add command -label "Use FFT" -command helpFFT
.f.help.m add command -label "How to learn" -command helpLearn

menubutton .f.fft -text "ft mode" -menu .f.fft.m -relief raised \
  -fg #ffffff -bg #c61e6e
menu .f.fft.m -fg #ffffff -bg #c61e6e
.f.fft.m add command -label "Help" -command helpFFT
.f.fft.m add command -label "Install FFT" -command fftInstall
.f.fft.m add command -label "Use FFT" -command { set UseFFT 1 }
.f.fft.m add command -label "Do Not Use FFT" -command { set UseFFT 0 }

button .f.smooth -text "smooth" -fg #ffffff -bg #b25f00 -command {
  if { $Smooth == 1 } {
    .f.smooth configure -text "smooth"
  } else {
    .f.smooth configure -text "polyln"
  }
  toggleSmooth
}

button .f.quit -text "quit" -command {exit} -fg #ffffff -bg #002332
#------------------------------------------------------------------------------
frame .fc
label .fc.mt -relief raised -width 40 \
  -fg #ffe6df -bg #8a0314 -text \
"Signal Waveform x(t)"

button .fc.pnt -text "hide pts" -fg #ffffff -bg #17714b -command {
  if { $ShowPnts == 1 } {
    .fc.pnt configure -text "show pts"
  } else {
    .fc.pnt configure -text "hide pts"
  }
  toggleShowPnts
}

### def func button ###
menubutton .fc.user -text "def func" -menu .fc.user.m -relief raised \
  -fg #ffffff -bg #6c5f44
menu .fc.user.m -fg #ffffff -bg #c61e6e
.fc.user.m add command -label "help" -command { helpUserfunc }
.fc.user.m add command -label "define a signal function" -command { userFunc }

### load data button ###
menubutton .fc.loaddata -text "load data" -menu .fc.loaddata.m -relief raised \
  -fg #ffffff -bg #a30a26
menu .fc.loaddata.m -fg #ffffff -bg #c61e6e
.fc.loaddata.m add command -label "help" -command { helpLoadsig }
.fc.loaddata.m add command -label "load a signal data" \
  -command { loadData }

### save data button ###
menubutton .fc.savedata -text "save data" -menu .fc.savedata.m -relief raised \
  -fg #ffffff -bg #a30a26
menu .fc.savedata.m -fg #ffffff -bg #c61e6e
.fc.savedata.m add command -label "help" -command { helpSavesig }
.fc.savedata.m add command -label "save the signal data" \
  -command { saveData }

### save eps button ###
menubutton .fc.epssig -text "save eps" -menu .fc.epssig.m -relief raised \
  -fg #ffffff -bg #b96e00
menu .fc.epssig.m -fg #ffffff -bg #c61e6e
.fc.epssig.m add command -label "help" -command { helpEpssig }
.fc.epssig.m add command -label "save the waveform in an eps file" \
  -command { epsSignal }

button .fc.clear -text "clear" -fg #ffffff -bg #5f1280 \
  -command {
  global N SpecIm SpecRe SpecPw Scale XYSpecIm XYSpecRe XYSpecPw \
    Resol X0 Y0 CurvIm CurvRe CurvPw
  for { set i 0 } { $i < $N } { incr i } {
    set Signal($i) 0; set SpecRe($i) 0; set SpecIm($i) 0; set SpecPw($i) 0
  }
  eval .c delete [.c find withtag line]
  eval .c delete [.c find withtag point]
  eval .d delete [.d find withtag line]
  ### draw signal waveform and plot control points ###
  set XYSignal [array2Coord Signal $N $Resol 1.0 $X0 $Y0]
  set CurvSig [eval .c create line $XYSignal -tags line -fill black]
  drawPoints .c Signal $N $Resol 1.0 $X0 $Y0
  ### draw spectrum (Re and Im) curves ###
  set Smooth 0
  set XYSpecIm [array2Coord SpecIm $N $Resol 1.0 $X0 $Y0]
  set CurvIm [eval .d create line $XYSpecIm -tags line -fill cyan -width 2]
  set XYSpecRe [array2Coord SpecRe $N $Resol 1.0 $X0 $Y0]
  set CurvRe [eval .d create line $XYSpecRe -tags line -fill violet -width 2]
  if { $ShowAbs == 1 } {
    set XYSpecPw [array2Coord SpecPw $N $Resol 1.0 $X0 $Y0]
    set CurvPw [eval .d create line $XYSpecPw -tags line -fill yellow]
  }
}
#------------------------------------------------------------------------------
frame .fd
label .fd.mf -relief raised \
  -fg #dfffff -bg #0d0564 -text \
"Fourier Spectrum X(f) :: violet: Re X(f), cyan: Im X(f), yellow: | X(f) |"
button .fd.abs -text "show abs" -fg #ffffff -bg #735500 -command {
  if { $ShowAbs == 1 } {
    .fd.abs configure -text "show abs"
  } else {
    .fd.abs configure -text "hide abs"
  }
  toggleShowAbs
}

### def func button ###
menubutton .fd.savespec -text "savedata" -menu .fd.savespec.m -relief raised \
  -fg #ffffff -bg #e64900
menu .fd.savespec.m -fg #ffffff -bg #c61e6e
.fd.savespec.m add command -label "help" -command { helpSavespec }
.fd.savespec.m add command -label "define a signal function" \
  -command { saveSpec }

### def func button ###
menubutton .fd.epsspec -text "save eps" -menu .fd.epsspec.m -relief raised \
  -fg #ffffff -bg #b96e00
menu .fd.epsspec.m -fg #ffffff -bg #c61e6e
.fd.epsspec.m add command -label "help" -command { helpSavespec }
.fd.epsspec.m add command -label "save the spectrum in an eps file" \
  -command { epsSpec }

button .fd.autoscale -text "man scale" \
  -command { toggleAutoScale } -fg #ffffff -bg #622662
#------------------------------------------------------------------------------
### scales ###
frame .s

# grid delta
scale .s.delta -label "time scale" -from 10 -to 400 \
  -length [expr ($W+13)/2] \
  -orient horizontal -command getDelta \
  -fg #ffffff -bg #355a41
pack .s.delta -side left -expand yes
proc getDelta value {
  global Delta Xorig W H TimeFreq
  set DivT [expr [.s.delta get]]
  eval .c delete [.c find withtag grid]

#  drawGrid $DivT
  for { set x $DivT } { $x < $W } { incr x $DivT } {
    set xx [expr $Xorig+$x]
    set id [.c create line $xx 0 $xx $H -fill grey50 -tags grid]
    .c lower $id
    set xx [expr $Xorig-$x]
    set id [.c create line $xx 0 $xx $H -fill grey50 -tags grid]
    .c lower $id
  }
  eval .d delete [.d find withtag grid]
#  empirically, 4600.0 was chosen to fit the relation: time x frequency = 1
#  set delta [expr 4600./$DivT]
  set delta [expr $TimeFreq/$DivT]
  for { set x $delta } { $x < $W } { set x [expr $x+$delta] } {
    set xx [expr round($Xorig+$x)]
    set id [.d create line $xx 0 $xx $H -fill grey50 -tags grid]
    .d lower $id
    set xx [expr round($Xorig-$x)]
    set id [.d create line $xx 0 $xx $H -fill grey50 -tags grid]
    .d lower $id
  }
  set GridShown 1
}
.s.delta set $DivT

# Spectrum display scaling factor (multipied 1/Scale)
scale .s.scale -label "spectrum scaling factor" -from 1 -to 50 \
  -length [expr ($W+13)/2] \
  -orient horizontal -command getScale \
  -fg #ffffff -bg #7d4930
pack .s.scale -side left -expand yes
proc getScale value {
  global N SpecIm SpecRe SpecPw Scale XYSpecIm XYSpecRe XYSpecPw ShowAbs \
    Resol X0 Y0 CurvIm CurvRe CurvPw
  set Scale [.s.scale get]
  set XYSpecIm [array2Coord SpecIm $N $Resol [expr 1.0/$Scale] $X0 $Y0]
  eval .d coord $CurvIm $XYSpecIm
  set XYSpecRe [array2Coord SpecRe $N $Resol [expr 1.0/$Scale] $X0 $Y0]
  eval .d coord $CurvRe $XYSpecRe
  if { $ShowAbs == 1 } {
    set XYSpecPw [array2Coord SpecPw $N $Resol [expr 1.0/$Scale] $X0 $Y0]
    eval .d coord $CurvPw $XYSpecPw
  }
}
.s.scale set $Scale
#------------------------------------------------------------------------------
#  packing and binding
#------------------------------------------------------------------------------
pack .fc.mt .fc.clear .fc.pnt .fc.user .fc.loaddata .fc.savedata .fc.epssig \
  -side left -fill both -expand yes
pack .fd.mf .fd.abs .fd.autoscale .fd.savespec .fd.epsspec \
  -side left -fill both -expand yes
pack .f.m .f.smooth .f.fft .f.help .f.quit \
  -side left -fill both -expand yes
pack .f .fd .d .fc .c .s -side top -fill both -expand yes

# dangerous - for quick debugging only
# bind . <Control-c> {destroy .}
# bind . <q> {destroy .}

focus .
#------------------------------------------------------------------------------
# point a point
bind .c <Button-1> {
  global This ThisY N ShowN
  set tags [lindex [.c itemconfigure current -tags] 4]
  if { [lsearch $tags point] >= 0 && [lsearch $tags current] >=0 } {
    .c itemconfigure current -fill red
    set curr [.c find withtag current]
    if { $This != $curr } { .c itemconfigure $This -fill yellow }
    set This $curr
    set ThisY [expr [lindex [.c coords $This] 1]+2]
  } else {
    .c itemconfigure $This -fill yellow
    set This "none"
#    .fc.value config -text ""
    .fc.mt config -text "Signal Waveform x(t)"
  }
}

# drag a point along the y-axis
bind .c <B1-Motion> {
  global This N ShowN
  if { $This != "none" } {
    set yy [expr [lindex [.c coords $This] 1]+2]
    .c move $This 0 [expr %y-$yy]
#    .fc.value config -text "x(t)=[expr $Yorig-%y]"
    .fc.mt config -text "Signal Waveform x(t)=[expr $Yorig-%y]"
  }
}

# release the button - DFT calculation
bind .c <ButtonRelease-1> {
  global This ThisY Signal SpecRe SpecIm Scale N Resol \
    X0 Y0 CurvSig XYSignal CurvIm XYSpecIm CurvRe XYSpecRe CurvPw XYSpecPw
  if { $This != "none" } {
    set yy [expr [lindex [.c coords $This] 1]+2]
    set i [lindex [lindex [.c itemconfigure $This -tags] 4] 1]
    set dsig [expr $ThisY-$yy]
    set Signal($i) [expr $Signal($i)+$dsig]
#    .fc.value config -text "x(t)=$Signal($i)"
    .fc.mt config -text "Signal Waveform x(t)=$Signal($i)"
    set j [expr 2*$i+1]
    set XYSignal [lreplace $XYSignal $j $j [expr $Yorig-$Signal($i)]]
    eval .c coord $CurvSig $XYSignal
  ### DFT incremental modification ###
    incrDFT $dsig $i SpecRe SpecIm SpecPw $N
    if { $AutoScale != 0 } {
      set Scale [expr [peakmax SpecRe $N]*2/$H]
      set scaleIm [expr [peakmax SpecIm $N]*2/$H]
      if { $scaleIm > $Scale } { set Scale $scaleIm }
      if { $ShowAbs == 1 } {
        set scaleIm [expr [peakmax SpecPw $N]*2/$H]
        if { $scaleIm > $Scale } { set Scale $scaleIm }
      }
    }
    set XYSpecIm [array2Coord SpecIm $N $Resol [expr 1.0/$Scale] $X0 $Y0]
    eval .d coord $CurvIm $XYSpecIm
    set XYSpecRe [array2Coord SpecRe $N $Resol [expr 1.0/$Scale] $X0 $Y0]
    eval .d coord $CurvRe $XYSpecRe
    if { $ShowAbs == 1 } {
      set XYSpecPw [array2Coord SpecPw $N $Resol [expr 1.0/$Scale] $X0 $Y0]
      eval .d coord $CurvPw $XYSpecPw
    }
  }
}

bind .c <ButtonRelease-2> { }
bind .c <ButtonRelease-3> { }

#------------------------------------------------------------------------------
# Fourier transform
#------------------------------------------------------------------------------
### discrete Fourier Transform (very slow) - not used here
proc DFT { signal specr speci n } {
  upvar $signal x
  upvar $specr zr
  upvar $speci zi
  set pi 3.14159265359
  for { set i 0 } { $i < $n } { incr i } {
    set sumr 0
    set sumi 0
    set domega [expr $pi*2*$i/$n]
    set omega 0
    for { set j 0 } { $j < $n } { incr j } {
      set omega [expr $omega+$domega]
      set sumr [expr $sumr+$x($j)*cos($omega)]
      set sumi [expr $sumi-$x($j)*sin($omega)]
    }
    set zr($i) [expr $sumr]
    set zi($i) [expr $sumi]
  }
}

### incremental DFT modification : when $signal($i) changed by $dsig
proc incrDFT { dsig i specr speci specp n } {
  upvar $specr zr
  upvar $speci zi
  upvar $specp pw
#  set pi 3.14159265359
#  set domega [expr $pi*2*($i-$n/2.)/$n]
  global OmegaMax
  set domega [expr $OmegaMax*($i-$n/2.)/$n]

#  set omega 0.0
  set omega [expr -$domega*$n/2]
  for { set j 0 } { $j < $n } { incr j } {
    set zr($j) [expr $zr($j)+$dsig*cos($omega)]
    set zi($j) [expr $zi($j)-$dsig*sin($omega)]
#    set pw($j) [expr sqrt($zr($j)*$zr($j)+$zi($j)*$zi($j))]
    set pw($j) [expr hypot($zr($j),$zi($j))]
    set omega [expr $omega+$domega]
  }
}

proc peakmax { array n } {
  upvar $array a
  set min -1
  set max 1
  for { set t 0 } { $t < $n } { incr t } {
    if { $a($t) > $max } { set max $a($t) }
    if { $a($t) < $min } { set min $a($t) }
  }
  set rev [expr -$min]
  if { $max < $rev } { return $rev }
  return $max
}
#------------------------------------------------------------------------------
# drawing
#------------------------------------------------------------------------------
### create array points - returns the point id list ###
proc drawPoints { wdw array n dx dy x y } {
  global Points
  upvar $array a
  for { set i 0 } { $i < $n } { incr i } {
    set xx [expr $i*$dx+$x]
    set yy [expr $y-$a($i)*$dy]
    set Points($i) \
      [$wdw create rectangle \
         [expr $xx-2] [expr $yy-2] [expr $xx+2] [expr $yy+2] \
         -fill yellow -tags [list point $i]]
  }
}

### convert an array into a coordinate list ###
proc array2Coord { array n dx dy x y } {
  upvar $array a
  set xy { }
  for { set i 0 } { $i < $n } { incr i } {
    lappend xy [expr $i*$dx+$x]
    lappend xy [expr $y-$a($i)*$dy]
  }
  return $xy
}

proc toggleSmooth { } {
  global Smooth CurvPw CurvIm CurvRe CurvSig ShowAbs
  if { $Smooth == 1 } {
    set Smooth 0
    .c itemconfig $CurvSig -smooth off
    if { $ShowAbs == 1 } { .d itemconfig $CurvPw -smooth off }
    .d itemconfig $CurvIm -smooth off
    .d itemconfig $CurvRe -smooth off
  } else {
    set Smooth 1
    .c itemconfig $CurvSig -smooth on
    if { $ShowAbs == 1 } { .d itemconfig $CurvPw -smooth on }
    .d itemconfig $CurvIm -smooth on
    .d itemconfig $CurvRe -smooth on
  }
}

proc toggleShowAbs { } {
  global SpecRe SpecIm SpecPw Scale N Resol X0 Y0 ShowAbs CurvPw XYSpecPw
  if { $ShowAbs == 1 } {
    set ShowAbs 0
    .d delete $CurvPw
  } else {
    set ShowAbs 1
#    powSpec SpecRe SpecIm SpecPw $N
    set XYSpecPw [array2Coord SpecPw $N $Resol [expr 1.0/$Scale] $X0 $Y0]
    set CurvPw [eval .d create line $XYSpecPw -tags line -fill yellow -width 1]
  }
}

proc toggleAutoScale { } {
  global AutoScale N H X0 Y0 Resol Scale ShowAbs \
    CurvRe CurvIm CurvPw XYSpecRe XYSpecIm XYSpecPw SpecRe SpecIm SpecPw
  if { $AutoScale } {
    set AutoScale 0
    .fd.autoscale configure -text "autoscale"
  } else {
    set AutoScale 1
    .fd.autoscale configure -text "man scale"
    if { $AutoScale != 0 } {
      set Scale [expr [peakmax SpecRe $N]*2/$H]
      set scaleIm [expr [peakmax SpecIm $N]*2/$H]
      if { $scaleIm > $Scale } { set Scale $scaleIm }
      if { $ShowAbs == 1 } {
        set scaleIm [expr [peakmax SpecPw $N]*2/$H]
        if { $scaleIm > $Scale } { set Scale $scaleIm }
      }
    }
    set XYSpecRe [array2Coord SpecRe $N $Resol [expr 1.0/$Scale] $X0 $Y0]
    eval .d coord $CurvRe $XYSpecRe
    set XYSpecIm [array2Coord SpecIm $N $Resol [expr 1.0/$Scale] $X0 $Y0]
    eval .d coord $CurvIm $XYSpecIm
    set XYSpecPw [array2Coord SpecPw $N $Resol [expr 1.0/$Scale] $X0 $Y0]
    eval .d coord $CurvPw $XYSpecPw
    update idletasks
  }
}

proc toggleShowPnts { } {
  global Signal N Resol X0 Y0 ShowPnts
  if { $ShowPnts == 1 } {
    set ShowPnts 0
    eval .c delete [.c find withtag point]
  } else {
    set ShowPnts 1
    drawPoints .c Signal $N $Resol 1.0 $X0 $Y0
  }
}
#------------------------------------------------------------------------------
# file I/O (save, load)
#------------------------------------------------------------------------------
proc epsSignal {} {
  toplevel .dialog
  label .dialog.label \
    -text "PostScript filename for Signal (click-3 to cancel):"
  entry .dialog.entry -relief sunken -textvariable fname
  pack .dialog.label .dialog.entry -side left -padx 1m -pady 1m
  bind .dialog.entry <Return> {
    .c postscript -colormode color -file $fname
    destroy .dialog
    grab set .
  }
  bind .dialog.label <ButtonPress-3> { destroy .dialog }
  bind .dialog.entry <ButtonPress-3> { destroy .dialog }
  grab set .dialog
}
proc epsSpec {} {
  toplevel .dialog
  label .dialog.label \
    -text "PostScript filename for Spectrum (click-3 to cancel):"
  entry .dialog.entry -relief sunken -textvariable fname
  pack .dialog.label .dialog.entry -side left -padx 1m -pady 1m
  bind .dialog.entry <Return> {
    .d postscript -colormode color -file $fname
    destroy .dialog
    grab set .
  }
  bind .dialog.label <ButtonPress-3> { destroy .dialog }
  bind .dialog.entry <ButtonPress-3> { destroy .dialog }
  grab set .dialog
}
proc saveData {} {
  global Signal
  toplevel .dialog
  label .dialog.label -text "Save signal data in (click-3 to cancel):"
  entry .dialog.entry -relief sunken -textvariable fname
  pack .dialog.label .dialog.entry -side left -padx 1m -pady 1m
  bind .dialog.entry <Return> {
    set file [open $fname w]
    for { set i 0 } { $i < $N } { incr i } {
      puts $file $Signal($i)
    }
    close $file
    destroy .dialog
    grab set .
  }
  bind .dialog.label <ButtonPress-3> { destroy .dialog }
  bind .dialog.entry <ButtonPress-3> { destroy .dialog }
  grab set .dialog
}
proc saveSpec {} {
  global Spec
  toplevel .dialog
  label .dialog.label -text "Save signal data in (click-3 to cancel):"
  entry .dialog.entry -relief sunken -textvariable fname
  pack .dialog.label .dialog.entry -side left -padx 1m -pady 1m
  bind .dialog.entry <Return> {
    set file [open $fname w]
    for { set i 0 } { $i < $N } { incr i } {
      puts $file $Spec($i)
    }
    close $file
    destroy .dialog
    grab set .
  }
  bind .dialog.label <ButtonPress-3> { destroy .dialog }
  bind .dialog.entry <ButtonPress-3> { destroy .dialog }
  grab set .dialog
}
proc loadData {} {
  global This ThisY Signal SpecRe SpecIm Scale N Resol ShowAbs \
    X0 Y0 CurvSig XYSignal CurvIm XYSpecIm CurvRe XYSpecRe CurvPw XYSpecPw
  toplevel .dialog
  label .dialog.label -text "Signal data file (click-3 to cancel):"
  entry .dialog.entry -relief sunken -textvariable fname
  pack .dialog.label .dialog.entry -side left -padx 1m -pady 1m
  bind .dialog.entry <Return> {
    set file [open $fname r]
    set i 0
    while { [gets $file newval] >= 0} {
#      .fc.value config -text "x(t)=$newval"
      .fc.mt config -text "Signal Waveform x(t)=$newval"
      set dsig [expr $newval-$Signal($i)]
      if { $dsig != 0.0 } {
        .c move $Points($i) 0 [expr -($dsig)]
        set Signal($i) $newval
        set j [expr 2*$i+1]
        set XYSignal [lreplace $XYSignal $j $j [expr $Yorig-$Signal($i)]]
        eval .c coord $CurvSig $XYSignal
        if { !$UseFFT } {
        ### DFT incremental modification ###
          incrDFT $dsig $i SpecRe SpecIm SpecPw $N
          if { $AutoScale != 0 } {
            set Scale [expr [peakmax SpecRe $N]*2/$H]
            set scaleIm [expr [peakmax SpecIm $N]*2/$H]
            if { $scaleIm > $Scale } { set Scale $scaleIm }
            if { $ShowAbs == 1 } {
              set scaleIm [expr [peakmax SpecPw $N]*2/$H]
              if { $scaleIm > $Scale } { set Scale $scaleIm }
            }
          }
          set XYSpecIm [array2Coord SpecIm $N $Resol [expr 1.0/$Scale] $X0 $Y0]
          eval .d coord $CurvIm $XYSpecIm
          set XYSpecRe [array2Coord SpecRe $N $Resol [expr 1.0/$Scale] $X0 $Y0]
          eval .d coord $CurvRe $XYSpecRe
          if { $ShowAbs == 1 } {
            set XYSpecPw [array2Coord SpecPw $N $Resol [expr 1.0/$Scale] $X0 $Y0]
            eval .d coord $CurvPw $XYSpecPw
          }
        }
        update idletasks
      }
      .c itemconfig $Points($i) -fill red
      incr i
    }
    close $file
    if { $UseFFT } {
      FFT Signal SpecRe SpecIm SpecPw $N
      if { $AutoScale != 0 } {
        set Scale [expr [peakmax SpecRe $N]*2/$H]
        set scaleIm [expr [peakmax SpecIm $N]*2/$H]
        if { $scaleIm > $Scale } { set Scale $scaleIm }
        if { $ShowAbs == 1 } {
          set scaleIm [expr [peakmax SpecPw $N]*2/$H]
          if { $scaleIm > $Scale } { set Scale $scaleIm }
        }
      }
      set XYSpecIm [array2Coord SpecIm $N $Resol [expr 1.0/$Scale] $X0 $Y0]
      eval .d coord $CurvIm $XYSpecIm
      set XYSpecRe [array2Coord SpecRe $N $Resol [expr 1.0/$Scale] $X0 $Y0]
      eval .d coord $CurvRe $XYSpecRe
      if { $ShowAbs == 1 } {
        set XYSpecPw [array2Coord SpecPw $N $Resol [expr 1.0/$Scale] $X0 $Y0]
      eval .d coord $CurvPw $XYSpecPw
      }
      update idletasks
    }
    destroy .dialog
    grab set .
    for { set i 0 } { $i < $N } { incr i } {
      .c itemconfig $Points($i) -fill yellow
    }
  }
  bind .dialog.label <ButtonPress-3> { destroy .dialog }
  bind .dialog.entry <ButtonPress-3> { destroy .dialog }
  grab set .dialog
}
#------------------------------------------------------------------------------
# user-defined function
#------------------------------------------------------------------------------
proc userFunc {} {
  global This ThisY Signal SpecRe SpecIm Scale N Resol ShowAbs \
    X0 Y0 CurvSig XYSignal CurvIm XYSpecIm CurvRe XYSpecRe CurvPw XYSpecPw
  toplevel .dialog
  label .dialog.label -text "Define a signal function (click-3 to cancel):"
  entry .dialog.entry -relief sunken -textvariable userfunc
  pack .dialog.label .dialog.entry -side left -padx 1m -pady 1m
  bind .dialog.entry <Return> {
    global This ThisY Signal SpecRe SpecIm Scale N Resol ShowAbs \
      X0 Y0 CurvSig XYSignal CurvIm XYSpecIm CurvRe XYSpecRe CurvPw XYSpecPw
#    .fc.user config -text "$userfunc"
    set funcshow [.c create text 10 10 -text "$userfunc" -anchor nw]
    for { set i 0 } { $i < $N } { incr i } {
      set t [expr ($i-$N/2)*$Resol/$DivT]
#      puts stdout [format "Time: %5.3f" $t]
      eval set newval [expr $userfunc]
#      .fc.value config -text "x(t)=$newval"
      .fc.mt config -text "Signal Waveform x(t)=$newval"
      update idletasks
      set dsig [expr $newval-$Signal($i)]
      if { $dsig != 0.0 } {
        .c move $Points($i) 0 [expr -($dsig)]
        set Signal($i) $newval
        set j [expr 2*$i+1]
        set XYSignal [lreplace $XYSignal $j $j [expr $Yorig-$Signal($i)]]
        eval .c coord $CurvSig $XYSignal
        if { !$UseFFT } {
        ### DFT incremental modification ###
          incrDFT $dsig $i SpecRe SpecIm SpecPw $N
          if { $AutoScale != 0 } {
            set Scale [expr [peakmax SpecRe $N]*2/$H]
            set scaleIm [expr [peakmax SpecIm $N]*2/$H]
            if { $scaleIm > $Scale } { set Scale $scaleIm }
            if { $ShowAbs == 1 } {
              set scaleIm [expr [peakmax SpecPw $N]*2/$H]
              if { $scaleIm > $Scale } { set Scale $scaleIm }
            }
          }
          set XYSpecIm [array2Coord SpecIm $N $Resol [expr 1.0/$Scale] $X0 $Y0]
          eval .d coord $CurvIm $XYSpecIm
          set XYSpecRe [array2Coord SpecRe $N $Resol [expr 1.0/$Scale] $X0 $Y0]
          eval .d coord $CurvRe $XYSpecRe
          if { $ShowAbs == 1 } {
            set XYSpecPw [array2Coord SpecPw $N $Resol [expr 1.0/$Scale] $X0 $Y0]
            eval .d coord $CurvPw $XYSpecPw
          }
        }
        update idletasks
      }
      .c itemconfig $Points($i) -fill red
    }
    if { $UseFFT } {
      FFT Signal SpecRe SpecIm SpecPw $N
      if { $AutoScale != 0 } {
        set Scale [expr [peakmax SpecRe $N]*2/$H]
        set scaleIm [expr [peakmax SpecIm $N]*2/$H]
        if { $scaleIm > $Scale } { set Scale $scaleIm }
        if { $ShowAbs == 1 } {
          set scaleIm [expr [peakmax SpecPw $N]*2/$H]
          if { $scaleIm > $Scale } { set Scale $scaleIm }
        }
      }
      set XYSpecIm [array2Coord SpecIm $N $Resol [expr 1.0/$Scale] $X0 $Y0]
      eval .d coord $CurvIm $XYSpecIm
      set XYSpecRe [array2Coord SpecRe $N $Resol [expr 1.0/$Scale] $X0 $Y0]
      eval .d coord $CurvRe $XYSpecRe
      if { $ShowAbs == 1 } {
        set XYSpecPw [array2Coord SpecPw $N $Resol [expr 1.0/$Scale] $X0 $Y0]
      eval .d coord $CurvPw $XYSpecPw
      }
      update idletasks
    }
    destroy .dialog
    .c delete $funcshow
    grab set .
    for { set i 0 } { $i < $N } { incr i } {
      .c itemconfig $Points($i) -fill yellow
    }
  }
  bind .dialog.label <ButtonPress-3> { destroy .dialog }
  bind .dialog.entry <ButtonPress-3> { destroy .dialog }
  grab set .dialog
}
#------------------------------------------------------------------------------
#  call external FFT algorithm
#------------------------------------------------------------------------------
proc FFT { signal specr speci specp n } {
  upvar $signal x
  upvar $specr zr
  upvar $speci zi
  upvar $specp zp
  if { $n!=128 } {
    puts stdout "!! Number of data ppoints must be 128 for use of FFF"
    exit
  }
 ### open "fft128" as an external function ###
  set p [open |fft128 r+]
  for { set i 0 } { $i < $n } { incr i } {
    puts $p $x($i)
  }
  flush $p
 ### receive the FFT results ###
  for { set i 0 } { $i < $n } { incr i } {
    gets $p zr($i)
  }
  for { set i 0 } { $i < $n } { incr i } {
    gets $p zi($i)
  }
  for { set i 0 } { $i < $n } { incr i } {
    gets $p zp($i)
  }
#  close $p
}
#------------------------------------------------------------------------------
# main program
#------------------------------------------------------------------------------
### draw axes ###
.c create line 0 $Yorig [expr $W+10] $Yorig -tags axis -fill grey20
.c create line $Xorig 0 $Xorig $H -tags axis -fill grey20
.d create line 0 $Yorig [expr $W+10] $Yorig -tags axis -fill grey80
.d create line $Xorig 0 $Xorig $H -tags axis -fill grey80

### draw signal waveform and plot control points ###
set XYSignal [array2Coord Signal $N $Resol 1.0 $X0 $Y0]
set CurvSig [eval .c create line $XYSignal -tags line -fill black -width 2]
drawPoints .c Signal $N $Resol 1.0 $X0 $Y0

### draw spectrum (Re and Im) curves ###
set XYSpecIm [array2Coord SpecIm $N $Resol 1.0 $X0 $Y0]
set CurvIm [eval .d create line $XYSpecIm -tags line -fill cyan -width 2]
set XYSpecRe [array2Coord SpecRe $N $Resol 1.0 $X0 $Y0]
set CurvRe [eval .d create line $XYSpecRe -tags line -fill violet -width 2]
#set XYSpecPw [array2Coord SpecPw $N $Resol 1.0 $X0 $Y0]
#set CurvPw [eval .d create line $XYSpecPw -tags line -fill yellow -width 1]

#------------------------------------------------------------------------------
# install an FFT program
#------------------------------------------------------------------------------
### write a C program and compile it ###
proc fftInstall {{w .basic}} {
  global FFTsrc
  set prog [open "/tmp/fft128.c" w]
  puts $prog $FFTsrc
  close $prog
  exec cc -s -o fft128 /tmp/fft128.c -lm
  exec rm /tmp/fft128.c
}
#------------------------------------------------------------------------------
set FFTsrc {\
/* fft128.c - an external program for "FourierTrans"
   coded by Shigeki Sagayama, 10 Oct 1998
   input: 128-point real data representing the center part of a 256-point data
   output: 128-point data with real and imaginary parts
*/

#include <stdio.h>
#include <math.h>
error(n) int n;
{
  fprintf(stderr,"?? number of input data is %d and not 128",n);
  exit(1);
}

main(argc,argv)
  int argc;
  char *argv[];
{
  int i,j;
  char line[80];
  double x[256],y[256],atof(),sqrt();

  for(i=0;i<256;i++) x[i]=y[i]=0; /* initialized */
  for(i=0;i<128;i++) {
    if(!gets(line)) error(i);
    x[(i+192)%256]=atof(line); /* the 64-th data will be stored in x[0] */
  }

  fft(x,y,8); /* 256-point FFT */

  for(i=0;i<128;i++) printf("%f\n",x[(i+192)%256]);
  for(i=0;i<128;i++) printf("%f\n",y[(i+192)%256]);
  for(i=0;i<128;i++) printf("%f\n",
    sqrt(x[(i+192)%256]*x[(i+192)%256]+y[(i+192)%256]*y[(i+192)%256]));
  fflush(stdout);
}

/* FFT program - original FORTRAN code in Oppenheim & Schafer: Digital Signal
   Processing, p. 331; corrected, converted into C and modified to speed it
   up by Shigeki Sagayama */
fft(xr,xi,m)
  double xr[],xi[];
  int m;
{
  int h,ih,i,j,k,l,n;
  double a,b,d,tr,ti,ur,ui,wr,wi;
  double sin(),cos(),pi=3.1415926535897932384626433832795;
  if(m<=0) return;
  n=1<<m; /* n=2^m; power */
/** radix-2 FFT **/
  if(m>=3) {
    h=n;
    for(l=m;l>2;l--) {
      k=h; h=k/2; ur=1.0; ui=0.0; d=pi/h;
      wr=cos(d); wi= -sin(d);
      for(j=0;j<h;j++) {
        for(i=j;i<n;i+=k) {
          ih=i+h; tr=xr[ih]; ti=xi[ih]; a=xr[i]-tr; b=xi[i]-ti;
          xr[ih]=a*ur-b*ui; xi[ih]=a*ui+b*ur; xr[i]+=tr; xi[i]+=ti; }
        a=ur*wr-ui*wi; ui=ur*wi+ui*wr; ur=a;
      }
  } }
/** 4-point DFT **/
  if(m>=2) {
    for(j=0;j<n;j+=4) {
      tr=xr[j]; ti=xi[j]; xr[j]=tr+xr[j+2]; xi[j]=ti+xi[j+2];
      xr[j+2]=tr-xr[j+2]; xi[j+2]=ti-xi[j+2];
      tr=xr[j+1]; ti=xi[j+1]; xr[j+1]=tr+xr[j+3]; xi[j+1]=ti+xi[j+3];
      tr=tr-xr[j+3]; ti=ti-xi[j+3]; xr[j+3]=ti; xi[j+3]= -tr;
  } }
/** 2-point DFT **/
  for(j=0;j<n;j+=2) {
    tr=xr[j]; xr[j]=tr+xr[j+1]; xr[j+1]=tr-xr[j+1];
    ti=xi[j]; xi[j]=ti+xi[j+1]; xi[j+1]=ti-xi[j+1];
  }
/** bit reversal **/
  h=n/2; j=0;
  for(i=0;i<n-1;i++) {
    if(i<j) {
      tr=xr[j]; ti=xi[j]; xr[j]=xr[i]; xi[j]=xi[i]; xr[i]=tr; xi[i]=ti; }
    k=h; while(k<=j) { j-=k; k/=2; } j+=k;
} }
}
#------------------------------------------------------------------------------
# end of a "wish" (tcl/tk) script: FourierTrans
#------------------------------------------------------------------------------