前回、scilab に搭載されている関数 FFT を使ってスペクトラム分析をしてみました。通常はデータ個数は 2 のべき乗個でなくてはいけない、とされているのですがなぜか scilab のはちゃんと動きます。本当に FFT 演算しているのか疑惑です。ただ処理速度を比較できる他のアプリケーションがないため、とりあえず信じるしかありません。
ただ今後はデータ個数は 2 のべき乗のみを扱うことにします。
さて、この FFT と今まで使っていた怪しげな手製スペクトラムアナライザとの実行時間の比較をしておくことにします。
実行時間の確認は関数 tic() toc() を使うのですが、マニュアルなどによりますとちょこちょこばらつくようです。そこで FFT の実行時間は同じことを 100回させてその平均時間を測定してみました。
スクリプトの全体はこんな感じです。データ個数は 2048 です。で、手製の方ですが本来なら半分の周波数までやれば良いはずなのですが、FFT が全周波数まで行っているため比較と云うことで同じ周波数までやっています。
xdel(winsid()) //ウィンドをクリア
clear
N=2048
t=0:1/N:1-1/N
p_shift=int(25/100*N) // 位相をずらす量をパーセントで記入
fig_axis=[0,1,-1.5,1.5]
Rect(1,1:N/4)=1
Rect(1,N/4+1:N/2)=1
Rect(1,N/2+1:N*3/4)=-1
Rect(1,N*3/4+1:N)=-1
if p_shift<>0 then
Srect(1:p_shift)=Rect(N-p_shift+1:N)
Srect(p_shift+1:N)=Rect(1:N-p_shift)
else
Srect=Rect
end
scf(1);clf;plot(t,Srect)
mtlb_axis(fig_axis)
// ここまでで波形を発生
for n=1:100
tic();
f_rect=fft(Srect) // 関数 fft を呼び出す。
time_duration(n)=toc()
M_rect=abs(f_rect)*2/N // 振幅を計算
end
disp(mean(time_duration))
// 手製スペクトラムアナライザ
tic();
fmax=N // 分析周波数範囲。実際はこの半分で良い。
for f1=1:fmax
x=sin(2*%pi*t*f1) // 正弦波
y=cos(2*%pi*t*f1) // 余弦波
combo_sin=x.*Srect
combo_cos=y.*Srect
sin_level(1,f1)=sum(combo_sin)./N*2
cos_level(1,f1)=sum(combo_cos)./N*2
level(1,f1)=((sin_level(f1))^2+(cos_level(f1))^2)^0.5
end
disp(toc())
fmax=20
scf(2);clf;plot(0:N/2-1,M_rect(1:N/2))
ft_axis=[0,fmax,0,1.5]
mtlb_axis(ft_axis)
結果は、disp で表示されているのですが、
FFT によるもの：0.9ms
手製：257ms
となりました。これでも数回やると多少変動します。
理論書などに書かれているものと直接比較する気はないですが、この結果でも分かるとおり FFT だけのことはあります。
ただ今後はデータ個数は 2 のべき乗のみを扱うことにします。
さて、この FFT と今まで使っていた怪しげな手製スペクトラムアナライザとの実行時間の比較をしておくことにします。
実行時間の確認は関数 tic() toc() を使うのですが、マニュアルなどによりますとちょこちょこばらつくようです。そこで FFT の実行時間は同じことを 100回させてその平均時間を測定してみました。
スクリプトの全体はこんな感じです。データ個数は 2048 です。で、手製の方ですが本来なら半分の周波数までやれば良いはずなのですが、FFT が全周波数まで行っているため比較と云うことで同じ周波数までやっています。
xdel(winsid()) //ウィンドをクリア
clear
N=2048
t=0:1/N:1-1/N
p_shift=int(25/100*N) // 位相をずらす量をパーセントで記入
fig_axis=[0,1,-1.5,1.5]
Rect(1,1:N/4)=1
Rect(1,N/4+1:N/2)=1
Rect(1,N/2+1:N*3/4)=-1
Rect(1,N*3/4+1:N)=-1
if p_shift<>0 then
Srect(1:p_shift)=Rect(N-p_shift+1:N)
Srect(p_shift+1:N)=Rect(1:N-p_shift)
else
Srect=Rect
end
scf(1);clf;plot(t,Srect)
mtlb_axis(fig_axis)
// ここまでで波形を発生
for n=1:100
tic();
f_rect=fft(Srect) // 関数 fft を呼び出す。
time_duration(n)=toc()
M_rect=abs(f_rect)*2/N // 振幅を計算
end
disp(mean(time_duration))
// 手製スペクトラムアナライザ
tic();
fmax=N // 分析周波数範囲。実際はこの半分で良い。
for f1=1:fmax
x=sin(2*%pi*t*f1) // 正弦波
y=cos(2*%pi*t*f1) // 余弦波
combo_sin=x.*Srect
combo_cos=y.*Srect
sin_level(1,f1)=sum(combo_sin)./N*2
cos_level(1,f1)=sum(combo_cos)./N*2
level(1,f1)=((sin_level(f1))^2+(cos_level(f1))^2)^0.5
end
disp(toc())
fmax=20
scf(2);clf;plot(0:N/2-1,M_rect(1:N/2))
ft_axis=[0,fmax,0,1.5]
mtlb_axis(ft_axis)
結果は、disp で表示されているのですが、
FFT によるもの：0.9ms
手製：257ms
となりました。これでも数回やると多少変動します。
理論書などに書かれているものと直接比較する気はないですが、この結果でも分かるとおり FFT だけのことはあります。