喫茶XRAY [Mr.XRAY掲示板]

[9115] Delphiで複素数計算

返信　削除

2015/8/24 (月) 13:04:37 taki

困っています。助けてください。

音声の線形予測分析でスペクトル包括を求めています。
ｎ次のLPC係数とそれぞれの残差分散が配列a[n],e[n]として求まった段階です。
スペクトル包括はFIRフィルター？のフィルタ係数として求めたa[n],e[n]を
代入すれば求まるそうです。

参考にしたWebページは
http://mitsui725.blogspot.jp/2013/08/c.html

C言語で書かれたプログラムをPascalに書き直しているのですが、複素数を扱う
complexというのが出てきます。
Delphi XE4 Starter ではVariant形式で複素数を扱うことができることが分かって
います。

悩んでいる内容は以下のような事です。

=== C言語 =========
double N2pi = 2*M_PI/N;
for (n=0;n<N + 1; ++n){
double complex z,x,・・略
x=-I*N2pi * n;
z=cexp(x);
　略

これを（我流で）書き換えると
var
  n:integer;
  N2pi,x,z:Variant;
begin
  N2pi := VarComplexCreate;
  x := VarComplexCreate;
  z := VarComplexCreate;
  N2pi := 2.0 * PI / NN;    // NNは分割する数で、外部からの値。
  for n := 0 to n do begin
    x := -I*N2pi * n;         // <- ここが分からない
    z := exp(x);

-I*N2pi * n;の意味を(-N2pi * n)を虚数部に代入する意味かなと思い、
x.Imaginary := -1.0 * N2pi * n; とかしたのですが、自身がありません。
実数部は0.0のままでいいのかな？

分かっている人からすれば、実際に数値を代入して計算させたら
正しいかどうか分かるんじゃない？と言われると思いますが、複素数の
計算をやったことが無いので、どんな計算をして確かめたらいいのかすら
分かりません。

何かヒントください。　m(_ _)m

環境は
Win7 Pro(32bit) + Delphi XE4 Starter

[9116] Re:Delphiで複素数計算

返信　削除

2015/8/24 (月) 14:22:19 福士

プログラミング言語 C の新機能
http://seclan.dll.jp/c99d/c99d09.htm

によると、C99のcomplexでは(大文字の)"I"は数学の虚数単位iのことだそうです。

あとStarter SKUに含まれているかどうかは判らないのですが、サンプルに
高度なレコード型を使ったcomplexの実装がありました。

C:\Users\Public\Documents\RAD Studio\11.0\Samples\Delphi\RTL\ComplexNumbers

あたりに見つかるとよいのですが…。参考にしてください。

[9117] Re2:Delphiで複素数計算

返信　削除

2015/8/24 (月) 16:54:26 taki

▼ 福士さん
> プログラミング言語 C の新機能
> http://seclan.dll.jp/c99d/c99d09.htm
>
> によると、C99のcomplexでは(大文字の)"I"は数学の虚数単位iのことだそうです。
>
> あとStarter SKUに含まれているかどうかは判らないのですが、サンプルに
> 高度なレコード型を使ったcomplexの実装がありました。
>
> C:\Users\Public\Documents\RAD Studio\11.0\Samples\Delphi\RTL\ComplexNumbers
>
> あたりに見つかるとよいのですが…。参考にしてください。

ありがとうございます。

Delphi付属のサンプルを探してみたのですが、見つけることができていません。

それと、z :=exp(x); は z := VarComplexExp(x); にしていました。

もう少し探してみます。付属のサンプルが確かにどこかにあった筈…。

[9118] Re3:Delphiで複素数計算

返信　削除

2015/8/24 (月) 17:24:29 taki

> もう少し探してみます。付属のサンプルが確かにどこかにあった筈…。

付属のサンプルはありました。
XEを起動してウェルカムページの「プロジェクト」より
「サンプルプロジェクトを開く」の中にありました。

しかし、Vassboth.Vcl.Compex.pasがあり、内容は
コンポーネントのようです。その他に〇〇.pasはありません。

Web上でも探してみます。

[9119] Re4:Delphiで複素数計算

返信　削除

2015/8/24 (月) 17:37:57 福士

> しかし、Vassboth.Vcl.Compex.pasがあり、内容は
> コンポーネントのようです。その他に〇〇.pasはありません。

これはコンポーネントではなく、普通のライブラリです。
高度なレコード型とオペレータオーバロードで複素数型を実装したものです。

[9121] Re5:Delphiで複素数計算

返信　削除

2015/8/24 (月) 19:40:02 taki

▼ 福士さん

> これはコンポーネントではなく、普通のライブラリです。
> 高度なレコード型とオペレータオーバロードで複素数型を実装したものです。

失礼しました。

じっくり見ると、色々な判別、関数計算がオーバーロードされていました。
IsZeroとかAbsとか、Expもありますね。

[9123] Re6:Delphiで複素数計算

返信　削除

2015/8/24 (月) 20:45:15 taki

procedure TForm1.Button3Click(Sender: TObject);
var
  a,b:Variant;
  rad,kaku:Double;    // 角度
begin
// 加算と減算が複素平面上でのベクトルの和と差で表され、
// 乗算と除算が複素平面上での回転

   a := VarComplexCreate;
   rad := (27 / 180) * pi;  // 27度
   a.real := cos(rad);      // 偏角27度を複素数にセット
   a.Imaginary := sin(rad);

   a := a * a;       // 複素数の積は偏角は和

   kaku := VarComplexAngle(a);   // radian
   kaku := kaku / pi * 180;
   label3.Caption := floattostr(kaku);   // 54 で正解

   // 偏角45度と偏角135(=180-45)度の足し算
   rad := (45 / 180) * pi;
   a.Real := cos(rad);
   a.Imaginary := sin(rad);

   rad := (135 / 180) * pi;
   b := VarComplexCreate(cos(rad),sin(rad));

   a := a + b;        // ベクトルの和

   kaku := VarComplexAngle(a);
   kaku := kaku / pi * 180;
   label4.Caption := floattostr(kaku);    // 90.000..1
   label5.Caption := floattostr(VarComplexAbs(a));   // 1.41421356...
end;

なんとなくわかってきました。

[9125] Re7:Delphiで複素数計算

返信　削除

2015/8/25 (火) 09:24:59 taki

C言語のコードの省略してある部分にヒントがありました。
denominator=numerator=0.0+I*0.0;

福士さんに教えてもらった"I"は数学の虚数単位iのことなので、
"+I*"より左の0.0が実数(Real)の値で、"+I*"より右の0.0が虚数
(Imaginary)の値のようです。

問題は、"+I*"じゃなくて"-I*"だということですね。

RealもImaginaryも実数を格納するから、やはり
0.0 +I* (-N2pi * n); の意味になるのかな？

今回参考にしたC言語のソースでは1.1+I*2.2の様式みたいですが、
1.1+2.2*Iのように記述する様式もあるみたいですね。

後は、カットアンドトライですか。道は険しい。

[9128] Re8:Delphiで複素数計算

返信　削除

2015/8/25 (火) 16:55:02 taki

// ------------------------------------

  rad := (27 / 180) * pi;  // 27度

  a.Real := rad;
  a.Imaginary := rad;

  kaku := VarComplexAngle(a);   // 偏角は？

  kaku := kaku / pi * 180;
  label1.Caption := floattostr(kaku);   // 常に結果は45

  label2.Caption := floattostr(VarComplexAbs(a));
  // rad := ( 27 / 180) * pi; の時は 0.6664...
  // rad := (127 / 180) * pi; の時は 3.134..   半径が変わる

// ------------------------------------

  rad := (127 / 180) * pi;

  a.Real := rad;            // 127
  a.Imaginary := rad;       // 127

  b := VarComplexExp(a);    // b = e^a

  kaku := VarComplexAngle(b);   // 偏角は?
  kaku := kaku / pi * 180;
  label3.Caption := floattostr(kaku);    // 127度
  label4.Caption := floattostr(VarComplexAbs(b));    // 9.1757...

// ------------------------------------

  rad := (127 / 180) * pi;  // 2.21656...

  a.Real := rad;            // 127
  a.Imaginary := -rad;      // -127

  b := VarComplexExp(a);    // b = e^a

  kaku := VarComplexAngle(b);   // 偏角は?
  kaku := kaku / pi * 180;
  label5.Caption := floattostr(kaku);    // -127度
  label6.Caption := floattostr(VarComplexAbs(b));     // 9.1757...

// ------------------------------------

色々と試してみました。

> 問題は、"+I*"じゃなくて"-I*"だということですね。

どうやら、右回りになりそうです。

> RealもImaginaryも実数を格納するから、やはり
> 0.0 +I* (-N2pi * n); の意味になるのかな？

Realが0なら、Imaginaryがどんな値をとっても・・・ですよね。

[9131] Re9:Delphiで複素数計算

返信　削除

2015/8/25 (火) 20:06:05 taki

// ------------------------------------

  rad := (127 / 180) * pi;

  a.Real := 0.0;            // 0.0
  a.Imaginary := rad;       // 127

  b := VarComplexExp(a);    // b = e^a

  kaku := VarComplexAngle(b);   // 偏角は?
  kaku := kaku / pi * 180;
  label3.Caption := floattostr(kaku);    // 127度
  label4.Caption := floattostr(VarComplexAbs(b));    // 1.0

// ------------------------------------

  rad := (127 / 180) * pi;  // 2.21656...

  a.Real := 0.0;            // 0.0
  a.Imaginary := -rad;      // -127

  b := VarComplexExp(a);    // b = e^a

  kaku := VarComplexAngle(b);   // 偏角は?
  kaku := kaku / pi * 180;
  label5.Caption := floattostr(kaku);    // -127度
  label6.Caption := floattostr(VarComplexAbs(b));     // 1.0

// ------------------------------------

> Realが0なら、Imaginaryがどんな値をとっても・・・ですよね。

すみません。Realが0でも意味がありました。
しかも、長さは1.0になりました。

[9132] そういうことか。

返信　削除

2015/8/25 (火) 23:08:18 taki

  rad := (127 / 180) * pi;

  a.Real := 0.0;            // 0.0
  a.Imaginary := rad;       // 127

  b := VarComplexExp(a);    // b = e^a

  kaku := VarComplexAngle(b);   // 偏角は?
  kaku := kaku / pi * 180;
  label3.Caption := floattostr(kaku);    // 127度
  label4.Caption := floattostr(VarComplexAbs(b));    // 1.0

// 改めて確認すると
  label5.Caption := floattostr(b.Real);          // -0.601815..
  label6.Caption := floattostr(b.Imaginary);     //  0.798635..

[9133] Re:そういうことか。

返信　削除

2015/8/26 (水) 19:02:47 Mr.XRAY

▼ takiさん

着々と進んでいるよですね．
なんにもお役に立てなくてスミマセン．

[9120] Re:Delphiで複素数計算

返信　削除

2015/8/24 (月) 18:01:44 Mr.XRAY

▼ takiさん，お疲れ様です．

> 困っています。助けてください。

助けられません (^^;

> 正しいかどうか分かるんじゃない？と言われると思いますが、複素数の
> 計算をやったことが無いので、どんな計算をして確かめたらいいのかすら
> 分かりません。

確かに．というより，数値計算の場合に言えますね．
昔，DEC のコンピュータ(パソコンじゃなくて)に付属していた数値計算
ライブラリをリファレンスにしていました．
そういうのがないと難しいのではないかと．

福士さんが書いている

>によると、C99のcomplexでは(大文字の)"I"は数学の虚数単位iのことだそうです。

これを見て，あっ! そうだ．c は大文字と小文字は区別するんだっけ．
と改めて認識しました．Delhpi 使っていると，区別ないですからね．
要注意です．

今回のとは直接関係があるかは分かりませんが．
数値計算関係が出てきたので，ついでに．
Delphi での情報は割と少ないようですが，皆さんの参考に．

[複素高速フーリエ変換]
http://www.w-frontier.com/delphi/tips.html#cfft

[岡本安晴のホームページ]
http://homepage2.nifty.com/okamotoy/

老舗ですが，これも
[TPMath] FFT のライブラリあり
https://www.unilim.fr/pages_perso/jean.debord/tpmath/tpmath.htm

どれも実務では使ったことがありません．あしからず．
この他にも数値計算関係のライブラリを見たことがあるのですが，失念．

INCM/CMT
Cyclamen v3.82