Octaveによるフーリエ解析

教科書の一連のフーリエ解析をoctaveで実施する。

octaveを使う利点
exp.d(波形データ)

#example wave ohsaki
#t d
0 5
0.5 32
1 38
1.5 -33
2 -19
2.5 -10
3 1
3.5 -8
4 -20
4.5 10
5 -1
5.5 4
6 11
6.5 -1
7 -7
7.5 -2

file:///home/shin/doc/lecture/hp/linux/expwave.png

(1)波形データの読み込み

> load exp.d  変数配列 exp に波形データが代入される。
> t=exp(:,1)  変数配列 t にexp の1列目のデータ(時刻)を代入する。
> d=exp(:,2)  
変数配列 d にexp の2列目のデータ(変位記録)を代入する。

(2)波形データの図示

> plot(t,d);

(3) フーリエ変換

> fd=fft(d)/16 フーリエ変換 16はデータの個数
> pd=fd.*conj(fd) パワースペクトルの計算 .*(要素同士のかけ算)とconj (共役複素数)
> tl= 0.5*16  時間長さ
> df=1/tl  周波数刻み
> f=(0:15)*df  周波数刻みの値をを行列で作成
> plot(f,abs(fd),"-;fd;",f,pd,"-;pd;")  パワースペクトル(pd)とフーリエ振幅(abs(fd))の図示
> print -landscape out.ps;  図をout.ps に出力する。

最後に ; (セミコロン)をつけると結果が画面に表示されない。



(4)平滑化

平滑化ではfor文を使ったマクロを組んだ。
平滑化のスクリプト hann.m

function ret=hann(x,nhan)
nhan
[n,m]=size(x)
ret=zeros(n,m);
temp=zeros(n,m);
nd=max(n,m);
F=[0.25,0.5,0.25];
temp=x;

for ihan=1:nhan
# ihan
ret(1)=0.5*(temp(1)+temp(2));
ret(nd)=0.25*temp(nd-1)+0.5*temp(nd)+0.25*temp(1);

for i=2:nd-1
ret(i)=F(1)*temp(i-1)+F(2)*temp(i)+F(3)*temp(i+1);
endfor
temp=ret;

endfor

スクリプトは,カレントフォルダーに置いておく。

> spd=hann(pd,3); hanninng windowsを3回かけ,結果をspdに代入する。
>  plot(f,spd,"+-;spd;",f,pd,"o-;pd;"); 結果の図示

(5)計算結果の保存

> out(:,1)=f;  fをoutの1列目に代入
> out(:,2)=pd; pdをoutの2列目に代入
> save out.d out;  outの値をout.dに保存する

(6)時刻歴での平滑化

> sd=hannt(d,2); hanning window で時刻歴波形を平滑化する

時刻歴波形を平滑化すると,高振動数領域をカットすることになる。

時刻歴波形を平滑化するスクリプト
function ret=hannt(x,nhan)
    nhan
    [n,m]=size(x)
    ret=zeros(n,m);
    temp=zeros(n,m);
    nd=max(n,m);
    F=[0.25,0.5,0.25];
        temp=x;
   
    for ihan=1:nhan
#          ihan
#          ret(1)=0.5*(temp(1)+temp(nd));
      ret(1)=0.25*temp(nd)+0.5*temp(1)+0.25*temp(2);
      ret(nd)=0.25*temp(nd-1)+0.5*temp(nd)+0.25*temp(1);

    for i=2:nd-1
          ret(i)=F(1)*temp(i-1)+F(2)*temp(i)+F(3)*temp(i+1);
    endfor
          temp=ret;
    endfor
    "hannt end"

平滑化された時刻歴波形と元波形のフーリエ振幅


(7)周波数応答関数を用いた応答計算

インパルス応答関数を計算するスクリプト
fresp.m
function ret=fresp(f,fr,h)

fr #固有振動数
h  #減衰定数
t1=f.*f;
t2=fr*fr-f.*f+i*2*h*f*fr;
ret=t1./t2;function ret=fresp(f,fr,h)

end

計算例
pulse.d

>load pulse.d
>t=pulse(:,1);
>d=pulse(:,2);
 
>load pulse.d
>t=pulse(:,1);
>d=pulse(:,2);
>dt=t(2)-t(1)
>size(t)
>nd=64
>tl=64*dt
>df=1/tl
>f=(0:63)*df
>f=f';
>fd=fft(d);
>tr=fresp(f,1.0,0.2)
>plot(f,tr)

>resp=tr.*fd
>resp=mkconj(resp)   折り返し振動数以上で共役複素数にする。

共役複素数を計算するスクリプト
mkconj.m
function ret=mkconj(x)

 [n,m]=size(x);
 nd=max(n,m)
 nfold=nd/2+1

 ret=zeros(n,m);
 temp=zeros(n,m);

 ret(1)=x(1);
 ret(nfold)=x(nfold);

for i=2:nfold-1
    k=2*nfold-i;
    ret(i)=x(i);
    ret(k)=conj(x(i));
endfor
endfunction


>respt=real(ifft(resp))
>plot(t,respt)