MATLABのfftでスペクトル解析

MATLABのfftでスペクトル解析

って言うほど大したことじゃないです。

fft使って周波数成分を抽出するだけです。

matlab_fft

けど、なれないとややこしいので紹介しますね。

ちなみに、上の図ですがMATLABのFFTサンプルから一部持ってきてwhite noizeってなってますが

duty比50%の矩形波なんですみません。

FFTとは

FFTとは高速フーリエ変換の略です。

連続系では単にフーリエ変換といいますが、コンピュータで処理する場合にデータを離散的に置きます。

それをDFTといい離散フーリエ変換の略ですが、DFTもコンピュータで処理できないわけじゃないと思いますが

アルゴリズムとして効率が悪いので、扱うデータを2のn乗のサンプルを用いて高速化するアルゴリズムです。

 

私の場合は、機械振動の周波数成分を抽出するのに使いますが、画像処理でもFFTってよく使うらしいです。

で、FFTは周波数解析に使うことがわかりました。

 

ではFFTのほかに周波数解析する手法なないのか?という話ですが・・・

私が知っているのはフィルター抽出です。

昔オーディオについていたグライコ(グラフィックイコライザー)や無線で使うディップメーターとか

そういった類のもので使います。

フィルターの場合は信号派に対して、フィルターの共振特性を使うものでして、

計測器の共振点をずらしながら共振周波数およびその大きさを取得するので時間がかかります。

が、よく使われます。

 

矩形波の周波数成分

矩形波は前回しめしたトランジスタアンプ的な発想で作りました。

周波数としては1kHzです。

図を見るとわかると思いますが、基本成分が1kHzをピークに

3kHz、5kHz・・・

とちょうど奇数倍のところにでてくるのがよくわかります。

 

矩形波というのは周期関数でして、周期関数というのはフーリエ級数で

示すことができるという特性があります。

 

そのフーリエ級数は式は省略しますが基本周波数の整数倍のサイン波形の

合成がフーリエ級数というものです。

 

FFTのスペクトル

グラフ見るとそのスペクトルは

1000Hz:1.273

3000Hz:0.4242

5000Hz:0.2543

7000Hz:0.1814

9000Hz:・・・

1100Hz:・・・

となります。

とりあえず1000~7000Hzまでのスペクトルを考えてグラフを描くとこんな感じです。

fft

 

w1は基本波のサイン波です

w2はw1に3倍成分乗せたものです

w3はw2に5倍成分乗せたものです

w4はw3に7倍成分乗せたものです

 

基本波に高調波を加えていくとどんどん矩形波に近くなってくことが

わかります。高調波を無限に考慮すれば完全な矩形波になることがわかります。

 

MATLABのプログラム

ここで使ったプログラムを紹介しますね。

SignalProcessingToolboxなどのMATLABオプションは不要です。

Fs=8192*2^4
aa=100*sin(2*pi*1000*[0:1/Fs:1]);
bb=2*gt(aa,1)-1;

subplot 211

plot([0:1/Fs:1],bb)
xlim([0 0.01])
grid on
title(‘squarewave’)
xlabel(‘time (milliseconds)’)
L=Fs;
NFFT = 2^nextpow2(L); % Next power of 2 from length of y
Y = fft(bb,NFFT)/L;
f = Fs/2*linspace(0,1,NFFT/2+1);

% Plot single-sided amplitude spectrum.
subplot 212
plot(f,2*abs(Y(1:NFFT/2+1)))
xlim([0 10000])
grid on
title(‘Single-Sided Amplitude Spectrum of y(t)’)
xlabel(‘Frequency (Hz)’)
ylabel(‘|Y(f)|’)

figure
tnn=[0:1/Fs:1];
w1=0.6366*2*sin(2*pi*1000*tnn);
w2=w1+0.4242*sin(2*pi*3*1000*tnn);
w3=w2+0.2543*sin(2*pi*5*1000*tnn);
w4=w3++0.1814*sin(2*pi*7*1000*tnn);
plot(tnn,w1,tnn,w2,tnn,w3,tnn,w4,’LineWidth’,2)
legend(‘w1′,’w2′,’w3′,’w4’)
hold on

xlim([0,0.002])
grid on

title(‘Fourier’)
xlabel(‘time (milliseconds)’)

sound(sin(2*pi*1000*[0:1/Fs:1]),Fs)
sound(bb,Fs)
sound(w4,Fs)

こんな感じですが、最後のsoundってのがmatlabのコマンドにあります。

引数のデータを音にしてスピーカーで鳴らしてくれます。

これもMATLABの特別なオプション使わないのでいいです。

 

Tags: ,

Leave a Reply