eof()は使うのに注意が必要(2)

例えば、あなたが4行からなるテキストファイルを持っていたとしましょう。

1 2.0
2 2.2
3 2.4
4 2.6

correct

// code correct
{ 
  ifstream finput("owondata3.txt");
  const Int_t icount_max=16;
  Int_t icount=0;
  Double_t data1[icount_max], data2[icount_max];
  
  while(icount < icount_max) { 
    finput >> data1[icount] >> data2[icount];
    std::cout << icount << ", " << data1[icount] << ", " << data2[icount] << std::endl;
    if(finput.eof()) break;
    icount++;
  }
  std::cout << "icount = " << icount << std::endl;

  TCanvas *c1 = new TCanvas("c1","Test",200,100,800,600);
  TGraph *g = new TGraph(icount, data1, data2);
  g->SetMarkerStyle( 20 );
  g->SetMarkerSize( 2.0 );
  g->Draw("AP");
}
% root correct.cc 
0, 1, 2
1, 2, 2.2
2, 3, 2.4
3, 4, 2.6
4, 0, 0
icount = 4
0, 1, 2
1, 2, 2.2
2, 3, 2.4
3, 4, 2.6

wrong

// code wrong
{
  ifstream finput("owondata3.txt");
  const Int_t icount_max=16;
  Int_t icount=0;
  Double_t data1[icount_max], data2[icount_max];

  while(! finput.eof() ) {
    finput >> data1[icount] >> data2[icount];
    std::cout << icount << ", " << data1[icount] << ", " << data2[icount] << std::endl;
    icount++;
  }
  std::cout << "icount = " << icount << std::endl;

  for(int i=0;i<icount;i++) {
    std::cout << i << ", " << data1[i] << ", " << data2[i] << std::endl;
  }

  TCanvas *c1 = new TCanvas("c1","Wrong",200,100,800,600);
  TGraph *g = new TGraph(icount, data1, data2);
  g->SetMarkerStyle( 20 );
  g->SetMarkerSize( 2.0 );
  g->Draw("AP");
}
% root wrong.cc 
0, 1, 2
1, 2, 2.2
2, 3, 2.4
3, 4, 2.6
4, 0, 0
icount = 5
0, 1, 2
1, 2, 2.2
2, 3, 2.4
3, 4, 2.6
4, 0, 0

eof()は使うのに注意が必要

もしも、あなたがAWKの1行プログラムを使わずに、ROOTのマクロで同じことをしようと思えば、コードは以下のようになります。

% cat owondata.txt | awk 'NR>8 {print $1, $3}' >owondata_ch1.txt
% cat owondata.txt | awk 'NR>8 {print $1, $4}' >owondata_ch2.txt
{
  ifstream finput("owondata.txt"); // read the original file without using AWK one-liners
  Int_t nskip=8;
  std::string s;

  for(int i=0;i<nskip;i++) { // skip the first eight lines
    getline(finput, s);
  }

  const Int_t icount_max=8192;
  Int_t icount=0;
  Double_t dummy, data1[icount_max], data2[icount_max], data3[icount_max]; // save $1, $3, and $4
  while(icount < icount_max) {
    finput >> data1[icount] >> dummy >> data2[icount] >> data3[icount];
    if(finput.eof()) break;
    icount++;
  }

  TCanvas *c1 = new TCanvas("c1","Test",200,100,800,600);
  TGraph *g = new TGraph(icount, data1, data2); // ch1 data
  g->SetMarkerStyle( 20 ); 
  g->SetMarkerSize( 0.5 );
  g->Draw("AP");
  TF1 *f1 = new TF1("f1", "[0]*sin([1]*x+[2])-[3]");
  f1->SetParameter(0, 2000.0);
  f1->SetParameter(1, 0.022);
  f1->SetParameter(2, 0.0);
  f1->SetParameter(3, 0.0);
  f1->SetLineColor(kRed);
  g->Fit(f1);

  TGraph *h = new TGraph(icount, data1, data3); // ch2 data
  h->SetMarkerStyle( 24 ); 
  h->SetMarkerSize( 0.5 );
  h->Draw("P");
  TF1 *f2 = new TF1("f2", "[0]*sin([1]*x+[2])-[3]");
  f2->SetParameter(0, 2000.0);
  f2->SetParameter(1, 0.022);
  f2->SetParameter(2, 0.0);
  f2->SetParameter(3, 0.0);
  f2->SetLineColor(kYellow);
  h->Fit(f2);

  double amp, phase;
  amp   = f2->GetParameter(0) / f1->GetParameter(0);
  phase = f2->GetParameter(2) - f1->GetParameter(2);
  std::cout << amp << ", " << phase << std::endl;

  TComplex v1, v2, v3, Z, z50;
  v1  = TComplex(1.0, 0.0);
  v2  = TComplex( amp*TMath::Cos(phase), amp*TMath::Sin(phase) );
  v3  = 2.0*v1 - v2;
  z   = v2/v3;
  z50 = 50.0*z;
  std::cout << v1 << ", " << v2 << ", " << v3 << ", " << z << ", " << z50 << std::endl;

  TComplex rho;
  double vswr;
  rho  = (z-1.0)/(z+1.0);
  vswr = (1.0+TComplex::Abs(rho)) / (1.0-TComplex::Abs(rho));
  std::cout << rho << ", " << vswr << std::endl;
}

以下のようなコードは、icountをwhile loopの外で1だけ減らさない限り、うまく動作しないことに留意して下さい。

while(! finput.eof() ) }
  finput >> ...
  icount++;
}

インピーダンス測定とカーブフィッティング(5)

従って、もし全てを1つのマクロで記述するとすれば、以下のようになります。

// file name = myetst41.cc
{
  TGraph *g = new TGraph("myfile3.dat");
  g->SetMarkerStyle( 20 ); 
  g->SetMarkerSize( 0.5 );
  g->Draw("AP");
  TF1 *f1 = new TF1("f1", "[0]*sin([1]*x+[2])-[3]");
  f1->SetParameter(0, 2000.0);
  f1->SetParameter(1, 0.022);
  f1->SetParameter(2, 0.0);
  f1->SetParameter(3, 0.0);
  f1->SetLineColor(kRed);
  g->Fit(f1);

  TGraph *h = new TGraph("myfile4.dat");
  h->SetMarkerStyle( 24 ); 
  h->SetMarkerSize( 0.5 );
  h->Draw("P");
  TF1 *f2 = new TF1("f2", "[0]*sin([1]*x+[2])-[3]");
  f2->SetParameter(0, 2000.0);
  f2->SetParameter(1, 0.022);
  f2->SetParameter(2, 0.0);
  f2->SetParameter(3, 0.0);
  f2->SetLineColor(kYellow);
  h->Fit(f2);

  double amp, phase;
  amp   = f2->GetParameter(0) / f1->GetParameter(0);
  phase = f2->GetParameter(2) - f1->GetParameter(2);
  std::cout << amp << ", " << phase << std::endl;

  TComplex v1, v2, v3, Z, z50;
  v1  = TComplex(1.0, 0.0);
  v2  = TComplex( amp*TMath::Cos(phase), amp*TMath::Sin(phase) );
  v3  = 2.0*v1 - v2;
  z   = v2/v3;
  z50 = 50.0*z;
  std::cout << v1 << ", " << v2 << ", " << v3 << ", " << z << ", " << z50 << std::endl;

  TComplex rho;
  double vswr;
  rho  = (z-1.0)/(z+1.0);
  vswr = (1.0+TComplex::Abs(rho)) / (1.0-TComplex::Abs(rho));
  std::cout << rho << ", " << vswr << std::endl;
}
% root mytest41.cc
0.778415, -0.53185
(1,0i), (0.670893,-0.394757i), (1.32911,0.394757i), (0.382788,-0.410701i), (19.1394,-20.535i)
(-0.329107,-0.394757i), 3.1148

ROOT、インタープリテーションかコンパイレーションか

ROOTには、ClingというC++インタープリターが内蔵されています。Clingを用いた1つの例は:

root6

% root
root [0] TF1 f1("f1", "sin(x)/x", -10.0, 10.0)
root [1] f1.Draw()

各行の最後に終端の”;”が不要なことに注意して下さい。

もしくは、あなたはROOTのmacro、これは行の集合ですが、を書くことができます。

// file name = myprog1.C
void myprog1(){
   TCanvas *c1 = new TCanvas("c1" ,"Test" ,0 ,0 ,600 ,400) ;
   c1->SetGrid( ) ;
   TF1 *f1 = new TF1("f1", "sin(x)/x", -10.0, 10.0);
   f1->Draw();
   c1->Update();
}

あなたは、あなたのマクロを以下のいずれかの方法で実行することができます。

% root myprog1.C
// or
% root
root [0] .x myprog1.C

root7

あなたのマクロは、いくらかの”dressing code”を付加することにより、コンパイルすることもできます。

https://root.cern.ch/root/htmldoc/guides/primer/ROOTPrimer.html#interpretation-and-compilation

// file name = myprog2.C

//include some header files, not necessary for Cling.
#include "TApplication.h"
#include "TROOT.h"
#include "TCanvas.h"
#include "TF1.h"
using namespace std;

void myprog2() {
  TCanvas *c1 = new TCanvas("c1", "Test", 0, 0, 600, 400);
  c1->SetGrid();
  TF1 *f1 = new TF1("f1", "sin(x)/x", -10.0, 10.0);
  f1->Draw();
  c1->Update();
}

void StandaloneApplication(int argc, char** argv) {
  myprog2();
}

int main (int argc, char** argv) {
  TApplication app("ROOT Application", &argc, argv);
  StandaloneApplication(app.Argc(), app.Argv());
  app.Run();
  return 0;
}
% g++ -o myprog2 myprog2.C `root-config --cflags --libs`
./myprog2

最後に、あなたはバインディングの集合PyROOTを用いてPythonにインターフェースすることもできます。

# file name = myprog3.py
from ROOT import TCanvas, TF1
c1=TCanvas("c1" ,"Test" ,0 ,0 ,600 ,400) 
c1.SetGrid ()
f1=TF1("f1", "sin(x)/x", -10.0, 10.0)
f1.Draw()
c1.Update()
raw_input('Press <ret> to end -> ')
% python myprog3.py

インピーダンス測定とカーブフィッティング(4)

vectors2

先ほど求めたパラメータp0とp2を用いれば、数値的には以下のようになります。

vector(red)    = 1.0*cexp(i*0.0) = (1.0, 0.0)
vector(green)  = vector(red) = (1.0, 0.0)
vector(yellow) = 0.778415*vector(red)*cexp(i*-0.531850) = (0.670893, -0.394757)
vector(purple) = 2*vector(red)-vector(yellow) = (1.32911, 0.394757)

これで、DUTのインピーダンスを計算する準備ができました。

Z(DUT) = 50 ohm * vector(yellow) / vector(purple)   <-- because the current is the same.
       = 50 ohm * (0.670893, -0.394757) / (1.32911, 0.394757)
       = 50 ohm * (0.382787,-0.4107)
       = (19.1394,-20.535) ohm

VSWRを計算するための短いプログラムです。

#include <complex>
#include <iostream>
using namespace std;
int main() {
 complex<double> x(0.670893, -0.394757);
 complex<double> y(1.32911 ,  0.394757);
 complex<double> rho;
 double vswr;
 cout << x << y << x/y << 50.0*x/y << endl;
 rho  = (x/y-1.0)/(x/y+1.0);
 vswr = (1.0+abs(rho))/(1.0-abs(rho));
 cout << rho << vswr << endl;
 return 0;
}

このプログラムの出力は:

% ./a.out
(0.670893,-0.394757)(1.32911,0.394757)(0.382787,-0.4107)(19.1394,-20.535)
(-0.329108,-0.394756)3.1148

Z(DUT)からVSWRを求めるのには、スミスチャートを用いることもできます。

smith1

インピーダンス測定とカーブフィッティング(3)

root4

vectors

図では、2つの信号がベクトルとして表現されています。赤色の方は基準信号で、黄色の方はDUTに掛かる信号です。実際には、2つのベクトルは信号源の周波数7026kHzで反時計回り(そう決めたのです)に回転していますが、私たちはベクトルの相対的な関係にしか興味がありません。

さて、どのようにしてDUTのインピーダンスを知ることができるでしょうか?測定回路のことを思い出せば、私たちは図中にさらに2つのベクトルを加えることができます。

directionalbridge2

vectors2

ここで、赤色と緑色のベクトルは同一であることに注意して下さい。なぜならば、2つの50オーム抵抗は直列に接続されているからです。紫色のベクトルは、赤色のベクトルを2倍したものと、黄色のベクトルとの差分として求めることができます。

インピーダンス測定とカーブフィッティング(2)

root5

さて、私のダイポールアンテナのインピーダンスを測定してみましょう。この図では、基準信号が赤色で、アンテナに加わる信号が黄色で表示されています。もしも、ダイポールアンテナのインピーダンスが正確に50オームであれば、2つの信号は同じになるはずです。

% head -15 owondata.txt
Save Time: 2016-09-11 13:27:26
Units:(mV)
                           CH1            CH2
Frequency:           7.026 MHz      7.042 MHz
Period:               0.142 uS       0.142 uS
SP:                   0.001 uS       0.001 uS
PK-PK:                 4.080 V        3.200 V

1              0.000                  1760.00        1560.00
2              0.001                  1760.00        1560.00
3              0.001                  1720.00        1560.00
4              0.002                  1720.00        1560.00

オシロスコープ、OWON PDS 5022S、から得られたテキストファイルは、ROOTのプログラムで処理される前に、AWKの1行プログラムで加工されます。

% cat owondata.txt | awk 'NR>8 {print $1, $3}' >owondata_ch1.txt
% cat owondata.txt | awk 'NR>8 {print $1, $4}' >owondata_ch2.txt
// file name = mytest5.cc
{
  TGraph *g = new TGraph("owondata_ch1.dat");
  g->SetMarkerStyle( 20 );
  g->SetMarkerSize( 0.5 );
  g->Draw("AP");
  TF1 *f1 = new TF1("f1", "[0]*sin([1]*x+[2])+[3]"); 
  f1->SetParameter(0, 2000.0);
  f1->SetParameter(1, 0.022);
  f1->SetParameter(2, 0.0);
  f1->SetParameter(3, 0.0);
  f1->SetLineColor(kRed);
  g->Fit(f1);

  TGraph *h = new TGraph("owondata_ch2.dat");
  h->SetMarkerStyle( 24 );
  h->SetMarkerSize( 0.5 );
  h->Draw("P");
  TF1 *f2 = new TF1("f2", "[0]*sin([1]*x+[2])+[3]");
  f2->SetParameter(0, 2000.0);
  f2->SetParameter(1, 0.022);
  f2->SetParameter(2, 0.0);
  f2->SetParameter(3, 0.0);
  f2->SetLineColor(kYellow);
  h->Fit(f2);
}

root4

% root mytest5.cc

****************************************
Minimizer is Minuit / Migrad
Chi2                      =       536853
NDf                       =          996
Edm                       =  7.63442e-09
NCalls                    =          148
p0                        =      2001.67   +/-   1.04511     
p1                        =    0.0220756   +/-   1.82267e-06 
p2                        =      2.03624   +/-   0.00102561  
p3                        =     -12.2567   +/-   0.749788    

****************************************
Minimizer is Minuit / Migrad
Chi2                      =       459226
NDf                       =          996
Edm                       =  8.91283e-08
NCalls                    =          128
p0                        =      1558.13   +/-   0.960465    
p1                        =    0.0220847   +/-   2.15852e-06 
p2                        =      1.50439   +/-   0.00124773  
p3                        =     -14.2113   +/-   0.68435     

パラメータp0とp2が、振幅と位相です。従って、振幅の比は1558.13/2001.67=0.778415、位相差は1.50439-2.03624=-0.531850[rad]=-30.4728[deg]です。

インピーダンス測定とカーブフィッティング

root3

あなたは被試験装置(DUT)、例えばあなたのアンテナ、のインピーダンスをそのデバイスに加えた正弦波の振幅と位相を観測することにより測定することができます。もちろん、測定は何らかの基準に対して行われるので、実際的にはあなたは2つの信号を同時に観測することになります。

owon7MHz20mOpen

例えばこの図では、あなたは赤色の信号を基準にして黄色の信号の振幅と位相を測定しようと思います。

詳細な議論については、私の頁Antennaを参照して下さい。

実際の測定は色々な方法で行うことができます。例えば、位相差を測定するために2つのゼロクロスポイント間でカーソルを手で動かすとか。ここでは、ROOTが備えるカーブフィッティング機能を用いた時にどうなるかを試してみようとしています。

短いC言語のプログラムを用いてテスト用のデータ・セットが生成されます。図では、黒色の丸と白色の丸とで表わされていて、赤色と黄色の曲線が得られた結果です。

% root -x mytest4.cc

****************************************
Minimizer is Minuit / Migrad
Chi2                      =  2.09673e-08
NDf                       =           47
Edm                       =   4.1926e-08
NCalls                    =          121
p0                        =      1.00001   +/-   4.08054e-06 
p1                        =     0.100001   +/-   3.45297e-07 
p2                        =      6.48313   +/-   8.9754e-06  

****************************************
Minimizer is Minuit / Migrad
Chi2                      =   9.1854e-10
NDf                       =           47
Edm                       =  1.83714e-09
NCalls                    =           70
p0                        =     0.500002   +/-   8.8495e-07  
p1                        =          0.1   +/-   1.29983e-07 
p2                        =      3.13999   +/-   3.40443e-06 
root [1] 

The parameters p0, p1, and p2 are the amplitude, the frequency, and the phase.

// file name = mytest4.cc
{
  TGraph *g = new TGraph("myfile4a.dat");
  g->SetMarkerStyle( 20 );
  g->SetMarkerSize( 1.0 );
  g->Draw("AP");
  TF1 *f1 = new TF1("f1", "[0]*sin([1]*x+[2])");
  f1->SetParameter(0, 0.25);
  f1->SetParameter(1, 0.11);
  f1->SetParameter(2, 3.14+0.2);
  g->Fit(f1);
  
  TGraph *h = new TGraph("myfile4b.dat");
  h->SetMarkerStyle( 24 );
  h->SetMarkerSize( 1.0 );
  h->Draw("P"); 
  TF1 *f2 = new TF1("f2", "[0]*sin([1]*x+[2])");
  f2->SetParameter(0, 0.25);
  f2->SetParameter(1, 0.11);
  f2->SetParameter(2, 3.14+0.2);
  h->Fit(f2);
}

バグキー(6)

root

ドット周期は時間の経過で変化するかと思ったのですが、何かを言うためにはもっと観察が必要です。

// file name = mytest.cc
{
  TCanvas *c = new TCanvas( "test" );
  TH1 *frame = c->DrawFrame( 0.0, 0.0, 11.0, 1200.0 );
  frame->GetXaxis()->SetTitle("time");
  frame->GetYaxis()->SetTitle("Duration");
  double x[11] = { 0.0, 1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0, 10.0 };
  double y[11] = { 1034, 883, 1028, 889, 1003, 923, 1053, 871, 1024, 901, 1048};
  TGraph *g = new TGraph( 11, x, y );
  g->SetMarkerStyle( 20 );
  g->SetMarkerColor(  4 );
  g->SetMarkerSize (  2 );
  g->Draw( "P" );
  TF1 *f1 = new TF1( "f1", "[0]*x+[1]" );
  gStyle -> SetOptFit(1111);
  g->Fit (f1);
}
% root -x mytest.cc
   -------------------------------------------------------------------------
  | Welcome to ROOT 6.06/08                             http://root.cern.ch |
  |                                            (c) 1995-2016, The ROOT Team |
  | Built for macosx64                                                      |
  | From heads/v6-06-00-patches@v6-06-06-30-g3bae07b, Sep 01 2016, 14:28:05 |
  | Try '.help', '.demo', '.license', '.credits', '.quit'/'.q'              |
   -------------------------------------------------------------------------

root [0] 
Processing mytest.cc...

****************************************
Minimizer is Minuit / Migrad
Chi2                      =      55135.1
NDf                       =            9
Edm                       =  2.85319e-19
NCalls                    =           32
p0                        =      1.30909   +/-   7.46271     
p1                        =      962.273   +/-   44.15       
root [1]