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

Programの内容

使用方法

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を利用していますので、任意のｘ軸に対してのフーリエ変換はできません。 ｘ軸についての情報はChange Wave Scaling を用いて設定してください。
3. このプログラムではデータ列の両端が0であるという前提で作られています。 一般のFFTでは演算を施す前に窓関数をかけたりします。また、位相のとびをなくすために 両端に0のデータ列を加えております。

/////////////////////////////////////////////////
// This procedure is for systematical fitting 
// of all waves in the active graph
//    created by M. Nagai on Jan., 2003,
//    modified on June. 30, 2009 
/////////////////////////////////////////////////

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
  $outf-=baseline    // set offset value
  InsertPoints 0,(pcsr(A)-1), $outf
  InsertPoints pcsr(B)+1,(2048-pcsr(B)), $outf
        
  FFT $outf 
  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)
        
  iterate(num)   // calculating power spectrum and phase
    mag=$outf[i]
    $outw[i]=cabs(mag)^2
    if (real(mag)<0)
        $outp[i]=atan(imag(mag)/real(mag))+Pi
    else
      $outp[i]=atan(imag(mag)/real(mag))
    endif
    if ($outp[i]<0)
      $outp[i]+=2*Pi
    endif
  loop
        
   // 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
        
  else
    abort "You inputed wrong wavename!"
  endif
endMacro