MATLABのfftでスペクトル解析
って言うほど大したことじゃないです。
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までのスペクトルを考えてグラフを描くとこんな感じです。
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の特別なオプション使わないのでいいです。