数値解析と制御のためのScilab入門

数値解析と制御のためのScilab入門では、行列、微分方程式、制御などの基礎的な計算を紹介しています。また、計算結果をScilabを用いたグラフィック表示による可視化について紹介しています。

*

高速フーリエ変換 FFTの計算

   

高速フーリエ変換 FFT の計算を行ってみます。

 フーリエ変換は雑音を含んだ N 個の離散データ からなる信号 x,x2,…,xn のなかからある周波数成分 Xを見つけるために用いられます。また、 X から x を求めることを逆変換といいます。

 X を効率よく計算するアルゴリズムが高速フーリエ変換(Fast Fourier Transform; FFT)と呼ばれています。scilabを用いて FFT 解析を行ってみることにします。

SPONSORED LINK

例として、サンプリング周波数を1kHzとし、周波数が50,90Hzの合成波のノイズがないものとして周波数を求めてみることにします。以下のようなプログラムを実行することによって得られます。


(1) ノイズなしの場合のFFT解析

1. // pure frequencies : 50 and 90 Hz
2. sample_rate=1000;
3. t = 0:1/sample_rate:0.5;
4. N=size(t,’*’); //number of samples(501)
5. clf();
6. s=sin(2*%pi*50*t)+sin(2*%pi*90*t+%pi/4);//+grand(1,N,’nor’,0,2);
7. subplot(2,1,1)
8. plot(t,s)
9. y=fft(s);
10. f=sample_rate*(0:(N/2))/N; //associated frequency vector
11. n=size(f,’*’)
12. subplot(2,1,2)
13. plot(f,abs(y(1:n)))


【プログラムの主なな概要】
2行目ではサンプリング周波数1000Hzとします。
3行目ではサンプリング時間を0.5(s)で時間ステップは0.001(s)です。
4行目ではサンプリング数で501個となります。
6行目でノイズがない場合の50,90Hzの合成波を求めています。
7行目でサブウインドウを指定しています。
8行目では合成波を表示しています。
9行目ではFFTの計算を行っています。
10行目で周波数を求めています。
11行ではサブウインドウを指定しています。
12行目ではfを横軸にFFTの値を表示します。


高速フーリエー変換 FFT

上の波形は、解析するための波形で、下の図はFFT解析による図です。
図からわかるように50Hz,と90Hzの波形のピーク値が見られます。周波数が混在した波形から、その波形が構成されている周波数を推測できます。
周波数解析においてFFTは、有効な手法の一つであることがわかります。


 - グラフィック, 組込み関数

        

SPONSORED LINK

SPONSORED LINK

  関連記事

カラーマップ:jetcolormap
3次元曲面をカラープロット surf で表示する

3次元グラフィックスを表示する場合に surf を使用します。また、プロットした …

ポリゴンを塗りつぶす
ポリゴンを任意の色で塗りつぶす xfpoly

ポリゴンを塗りつぶしたい場合があるかと思いっます。 塗りつぶす命令として xfp …

3次元グラフィックス meshgridとmesh

表示する場合には組み込み関数eshgrid, mesh をついで使用します。 そ …

2次元プロットにグリッドを追加
2次元プロットにグリッドを追加:xgrid|グラフィック

グラフのデータを比較したいとか、見やすくしたい場合には、グリッドを入れて表示させ …

グラフィック・ウインドウに
複数の色を用いてグラフを描く plot2d

ヒトツノウインドウにラインの色を変化させて、数個のグラフを重ねて表示することにし …

scilab 小数点以下切り捨て
小数点以下切り捨て floor|scilab入門

小数点以下を切り捨てる場合に:floor を使います。 例として  1.523 …

sinのグラフ
scilabのエディター SciNoteの使用

プログラムの作成に「SciNoteを起動」を使用するととても便利です。プログラム …

scilab 絶対値
絶対値の計算 abs|scilab入門

絶対値の計算:abs の計算の紹介です。  5 – 9 の計算は   …

scilab 切り上げ
小数点以下切り上げ ceil | scilab入門

小数点以下切り上げる場合には、ceil を満ちいます。 例えば  0.5 の場合 …

scilab パイi
円周率 π を表す|scilab入門

円周率 π を表す場合には、パーセントの後にpiとし  %pi のように表記しま …