データ解析のためのIgorPro活用術

電場波形からパワースペクトルと位相を計算する

このIgorPro のProcedureは、テラヘルツ時間領域分光法 などで得られた 電場波形からパワースペクトルと位相を計算するためのプログラムです。

Programの内容

このプログラムではwavenameという名のwaveをまずwavename_FFT という名前で複製し、データの前後に2000Points ほどのデータをInsertします。 そしてこのデータをFFTをかけます。得られたwavename_FFTは 複素数ですので、パワースペクトル(絶対値の2乗)をに、位相をとして計算させます。 このとき、位相は-PiからPiまでの範囲ですので、位相のとびが おこならないような演算を施します。

使用方法

  1. 演算させたいwaveをグラフ上に表示させ、そのグラフをActiveWindowにする。 また、演算させたい範囲をカーソルを用いて指定する。
  2. Procedure Windowにpowerspectrumを貼り付ける。なお、そのときに前後にInsertする Pointの数を修正してもよい。
  3. powerspectrum()を実行し、演算させたいWaveを選択する。
  4. 最終的にはwavename_FFTに複素振幅、wavename_Iにパワースペクトル、 wavename_pに位相のwaveが出力されます。

カスタマイズ、その他

  1. このプログラムのKeyとなるコマンドはもちろんFFTとなっております。 また、複素演算においてはcabs(z), real(z), imag(z)が重要な関数となります
  2. FFTを利用していますので、任意のx軸に対してのフーリエ変換はできません。 x軸についての情報はChange Wave Scaling を用いて設定してください。
  3. このプログラムではデータ列の両端が0であるという前提で作られています。 一般のFFTでは演算を施す前に窓関数をかけたりします。また、位相のとびをなくすために 両端に0のデータ列を加えております。

2022.6.23 普段使っているマクロに変えました。以前のマクロはこちらから
ここからProcedureです。

Macro powerspectrum(inw, baseline)
  String inw
  Variable/D baseline=0
  Prompt inw, "select wave",popup WaveList("*", ";","WIN:")
  Prompt baseline, "input baseline"
  String outf=inw+"_FFT"
  String outw=inw+"_I"
  String outp=inw+"_p"
  String outc=inw+"_c"
  Variable/C mag
  Variable num,num2
  silent 1
  PauseUpdate
        
  if (exists(inw)==1)
  Duplicate/R=(xcsr(A),xcsr(B))/O $inw $outf
  if(mod(numpnts($outf),2)==0)      
  FFT $outf
  else
  Duplicate/R=(xcsr(A),xcsr(B)+deltax($inw))/O $inw $outf  
  FFT $outf
  endif 
  num=numpnts($outf)
  Make/N=(num)/D/O $outw, $outp
  SetScale/P x 0,(1/num/deltax($inw)/2),"", $outw, $outp
  print xcsr(A), xcsr(B), deltax($inw)
        
  $outw=real(r2polar($outf))^2
  $outp=imag(r2polar($outf))
   Unwrap 2*pi, $outp
   
     // phase control
  Duplicate/O $outp $outc
  iterate(num-1)
    num2=num-1-i
    $outc[num2]=$outp[num2]-$outp[num2-1]
    if ($outc[num2]>4)
      $outc[num2]=0
      $outc[num2,num-1]-=2*Pi
      else
      if ($outc[num2]<-0.8)
        $outc[num2]= 0
        $outc[num2,num-1]+= 2*Pi
      else
        $outc[num2]=0
      endif
    endif
  loop
  
  $outp+=$outc
  Killwaves $outc
   
   $outp+=2*pi*x*xcsr(A)

  else
    abort "You inputed wrong wavename!"
  endif
endMacro


ここまで
最終更新日: 2004.4.1
一つ上の項目に戻る