円柱の周りの流れを解析しようとしており、風洞実験から得た一連の Cp 値があります。最初に、20 Hz のサンプル周波数から始めて、matlab で FFT を使用して渦放出の周波数を見つけようとしました。約7 Hzの周波数が得られました。次に、同じ実験を行いましたが、サンプリング周波数を 20 Hz から 200 Hz に変更しただけです。渦放出の周波数は約 70 Hz になりました (これがグラフのピーク位置です)。入力した Cp データに関係なく、グラフは変化しません。ピークが異なるのは、サンプル周波数を変更したときだけです。渦放出の頻度の増加はサンプル頻度に比例しているように見えますが、これはまったく意味がないようです。
7 に答える
あなたが見ている問題は、ナイキスト周波数(サンプリング周波数の半分)よりも高い周波数を検出できるFFTの制限による「データエイリアシング」に関連しています。
データエイリアシングを使用すると、実際の周波数のピークが中心になります(ナイキスト周波数を法とする実際の周波数)。20 Hzのサンプリングでは(70 Hzが真の周波数であると仮定すると、周波数がゼロになり、実際の情報が表示されないことを意味します。これに役立つ1つのことは、FFT「ウィンドウ処理」を使用することです。
発生している可能性のあるもう1つの問題は、単一FFT測定によるノイズの多いデータ生成に関連しています。大量のデータを取得し、オーバーラップのあるウィンドウ処理を使用し、結果を見つけるために平均するFFTが少なくとも5つあることを確認することをお勧めします。Steven Loweが述べたように、可能であれば、より速いレートでサンプリングする必要もあります。楽器がサンプリングできる最速のレートでサンプリングすることをお勧めします。
最後に、 Cの数値レシピからの抜粋を読むことをお勧めします(<-リンク):
- セクション12.0-FFTの概要
- セクション12.1(データエイリアシングについて説明します)
- セクション13.4(FFTウィンドウ処理について説明)
Cのソースコードを読む必要はありません。説明だけです。Cの数値レシピには、この主題に関する優れた要約情報があります。
他にご不明な点がございましたら、コメント欄にご記入ください。私は彼らに答えるために最善を尽くします。
幸運を!
これはおそらくプログラミングの問題ではなく、実験と測定の問題のように思えます
サンプリング周波数は、発振周波数の少なくとも 2 倍のレートである必要があると思います。そうしないと、アーティファクトが発生します。これは違いを説明するかもしれません。どちらの場合も、サンプリング周波数に対する FFT 周波数の比率は 0.35 であることに注意してください。より高いサンプリング レートで実験を繰り返すことができますか? これが強風の中で狭いシリンダーである場合、サンプリング レートが検出できるよりも速く振動/振動している可能性があると考えています..
これがお役に立てば幸いです - 97.6% の確率で、私が何について話しているのかわからない可能性があります ;-)
エイリアシングの問題でない場合は、正規化された周波数スケールで周波数応答をプロットしているように聞こえますが、これはサンプル周波数によって変化します。以下は、Matlab で信号の周波数応答をプロットする適切な方法の例です。
Fs = 100;
Tmax = 10;
time = 0:1/Fs:Tmax;
omega = 2*pi*10; % 10 Hz
signal = 10*sin(omega*time) + rand(1,Tmax*Fs+1);
Nfft = 2^8;
[Pxx,freq] = pwelch(signal,Nfft,[],[],Fs)
plot(freq,Pxx)
pwelch
「実際の」周波数データを出力するには、サンプル周波数をコマンドに明示的に渡す必要があることに注意してください。そうしないと、サンプル周波数を変更すると、共鳴が発生するビンがシフトするように見えます。これは、説明した問題に似ています。
DFT(FFT)のすべてのニュアンスを理解し始める前に、デジタル信号処理について真剣に読む必要があると思います。もし私があなただったら、まずこの素晴らしい本を読んで基礎を固めるだろう:
あなたの能力を本当に拡張する数学的処理がもっと必要な場合は、
この関連する質問を見てください。元々は VB について質問されていましたが、応答は一般的に FFT に関するものです
上記の周波数応答コードを使用してみましたが、Matlab に適切なツールボックスがないようです。fft コマンドを使用せずに同じことを行う方法はありますか? これまでのところ、これは私が持っているものです:
% FFT Algorithm
Fs = 200; % Sampling frequency
T = 1/Fs; % Sample time
L = 65536; % Length of signal
t = (0:L-1)*T; % Time vector
y = data1; % Your CP values go in this vector
NFFT = 2^nextpow2(L); % Next power of 2 from length of y
Y = fft(y,NFFT)/L;
f = Fs/2*linspace(0,1,NFFT/2);
% Plot single-sided amplitude spectrum.
loglog(f,2*abs(Y(1:NFFT/2)))
title(' y(t)')
xlabel('Frequency (Hz)')
ylabel('|Y(f)|')
使用しているコードに問題がある可能性があると思います。よくわかりませんが。
私の同僚は、スペクトル分析用のいくつかの素晴らしい GPL ライセンス関数を作成しました: http://www.mecheng.adelaide.edu.au/~pvl/octave/
(更新: このコードは現在、Octave モジュールの 1 つの一部です
:
http://octave.svn.sourceforge.net/viewvc/octave/trunk/octave-forge/main/signal/inst/。
そこから必要な部分だけを抽出します。)
これらは Matlab と Octave の両方向けに作成されており、主に Signal Processing Toolbox の類似関数のドロップイン置換として機能します。(したがって、上記のコードは引き続き正常に動作するはずです。)
データ分析に役立つ場合があります。などで自分で転がすよりも優れていますfft
。