パワースペクトルからの時系列信号生成

質問

独学でMatlabを用いた信号処理を学んでいます。
ある論文で、パワースペクトルを所与として時系列信号を生成するということをしています。再現したいのですが、その具体的な実装方法がわからず困っています。

以下に自分の理解を書きます。
いまある時系列（長さn;2の階乗）をy、そのfft結果をY、そしてパワースペクトルをPyと表します。通常の解析では、パワースペクトルは
　Py = Y.*conj(Y) / n
で求められると思います。いまは逆にPyが与えられています。Yからyに戻すのは逆フーリエ変換により可能です。しかし、PyからYへの戻し方がわかりません。単にPyのルートをとるのでは複素数情報が失われてしまい結果が違ってしまいます。
なにかとても簡単なところでつまづいているだけな気はするのですが、ご教示いただければ幸いです。どうぞよろしくお願いいたします。

参考までに、論文ではPyはたとえば次のような式です。
Py = 0.625 ./ ((abs( (exp(-j.*2.*pi.*f.*dt) - 1.02 .* exp(-20.*j.*pi.*dt)) ...
.* (exp(-j.*2.*pi.*f.*dt) - 1.02 .* exp(20.*j.*pi.*dt)) )).^2);

回答No.2

gukky · ベストアンサー

#1です。
前回の回答の式で一部間違えていました。
誤：y(t)=Σ(a(i)*cos(p(i)))
正：y(t)=Σ(a(i)*cos(ω*i*p(i)))
ωは対象信号の基本角周波数です。

>1．お書きいただいた「振幅情報」は、具体的にパワースペクトルからどのように求めればよいのでしょうか？
振幅情報＝パワースペクトラムの平方根です。

>2．論文には「位相情報」については一切ふれられていません。この場合、なにか暗黙の了解のようなものがあるのでしょうか。何もない場合、これだけでは論文の結果が再現できないことになってしまうのでしょうか…。
位相情報がない場合、元の時系列データは再現できません。

パワースペクトラムは、ある周波数成分がどれくらいの量含まれているかを示すものです。

簡単な例として下記の場合を考えます。
y=0.5*sin(ω*t+p1)+0.2*sin(2*ω*t+p2)+1.5*sin(3*ω*t+0.5*π+p3)
この時系列データでは基本周波数とその2倍及び3倍の高調波成分が含まれたものです。
これにfftを行うと、パワースペクトラムの値はp1,p2,p3の値が変わっても同じ値になります。
しかし、p1,p2,p3が変わると時系列データは変わります。
ですから、このp1,p2,p3である位相情報がないと元の時系列データを特定することはできません。

お礼

よくわかりました。
ご丁寧な解説、どうもありがとうございました。
お礼が遅れてしまったこと、お詫び申し上げます。

回答No.1

gukky

ある時系列データをfftすると、振幅情報と位相情報が得られます。パワースペクトルはその振幅情報のみから得られる結果なので、これだけでは元の波形を作ることはできません。
つまり、位相情報がないと駄目です。
matlabでの記述方法は覚えていないのですが、一般的には以下の式で表せばよいはずです。
y(t)=Σ(a(i)*cos(p(i)))
ここで、a(i)は振幅情報で、p(i)は位相情報です。
で、iは０～n/2です。
尚、通常fft解析結果の位相情報は余弦で与えるようですので、ここではcos(p(i))としています。

補足

gukkyさま、ありがとうございます。

申し訳ないのですが、正直まだ理解しきれておりません。以下の2点ご教示いただければ幸いです。自分がやりたいのは、与えられたパワースペクトルの式から模擬時系列を発生させることです。

1．お書きいただいた「振幅情報」は、具体的にパワースペクトルからどのように求めればよいのでしょうか？
　(Py * n)のルート を逆フーリエ変換
のでよいのでしょうか？
（計算すると複素数が出てきてしまいますが…）

2．論文には「位相情報」については一切ふれられていません。この場合、なにか暗黙の了解のようなものがあるのでしょうか。何もない場合、これだけでは論文の結果が再現できないことになってしまうのでしょうか…。