4

周期的なデータのセットがあります (正弦波ではありません)。1 つのベクトルに時間値のセットがあり、2 番目のベクトルに振幅のセットがあります。関数の周期をすばやく概算したいと思います。助言がありますか？

具体的には、これが私の現在のコードです。ベクトル t に対するベクトル x(:,2) の周期を概算したいと思います。最終的には、多くの初期条件に対してこれを行い、それぞれの周期を計算して結果をプロットしたいと思います。

function xdot = f (x,t)
         xdot(1) =x(2);
         xdot(2) =-sin(x(1));
endfunction

x0=[1;1.75];     #eventually, I'd like to try lots of values for x0(2)
t = linspace (0, 50, 200);


x = lsode ("f", x0, t)

plot(x(:,1),x(:,2));

ありがとうございました！

ジョン

4

2 に答える 2

6

自己相関関数を見てみましょう。

ウィキペディアより

自己相関は、信号とそれ自体の相互相関です。非公式には、観測間の時間間隔の関数としての観測間の類似性です。これは、ノイズに埋もれた周期信号の存在などの繰り返しパターンを見つけたり、その高調波周波数によって暗示された信号の欠落した基本周波数を特定したりするための数学的ツールです。これは、時間領域信号などの関数または一連の値を分析するための信号処理でよく使用されます。

Paul Bourke は、高速フーリエ変換に基づいて効率的に自己相関関数を計算する方法について説明しています (リンク)。

于 2010-04-04T08:43:36.010 に答える
2

離散フーリエ変換は、周期性を与えることができます。時間枠を長くすると周波数分解能が高くなるため、t定義をに変更しましたt = linspace(0, 500, 2000)。 時間領域 http://img402.imageshack.us/img402/8775/timedomain.png (ここにプロットへのリンクがあります。ホスティング サイトの方がよく見えます)。あなたがすることができます：

h = hann(length(x), 'periodic'); %# use a Hann window to reduce leakage
y = fft(x .* [h h]); %# window each time signal and calculate FFT
df = 1/t(end); %# if t is in seconds, df is in Hz
ym = abs(y(1:(length(y)/2), :)); %# we just want amplitude of 0..pi frequency components
semilogy(((1:length(ym))-1)*df, ym);

周波数ドメイン http://img406.imageshack.us/img406/2696/freqdomain.png プロット リンク。

グラフを見ると、最初のピークは約 0.06 Hz にあり、これは に見られる 16 秒の周期に対応しplot(t,x)ます。

ただし、これは計算上それほど高速ではありません。FFT は N*log(N) 操作です。

于 2010-04-04T08:18:53.413 に答える