前の記事、
FFT 使い手 Level 5 にありがちなこと（ランプ関数の影響を少なくするには）
http://blogs.yahoo.co.jp/susanoo2001_hero/12412381.html
では、手製のスペクトラムアナライザでランプ関数の影響を少なくするための実験を行ってみました。
今度は scilab の FFT で同じことをやってみます。
スクリプトはこれです。
FFT の区間は 1.0 ～ 2.0 秒です。 10Hz 以外は DC が乗ってしまいます。まあこれは結果の 0Hz を無視すると云うことでいいでしょう。
xdel(winsid()) //ウィンドをクリア
clear // Data をクリア
t_range=2 // 測定時間幅
N=2048 // サンプリング周波数
t=0:1/N:t_range-1/N // 時間内のサンプル数
fig_axis=[0,2,-2.0,2.0]
fs=10 // 正弦波の周波数
amp_sin=0.1 // 正弦波の振幅
Ramp(1:N*t_range)=[-(N*t_range)/2+0.5:1:(N*t_range)/2-0.5]/N*2 // ランプ関数
phi=100/N // 正弦波の位相
S_sin(1,1:N*t_range)=sin(2*%pi*(t+phi)*fs)*amp_sin // 本来の正弦波
Signal=Ramp+S_sin
scf(10);clf;plot(t,Signal)
mtlb_axis(fig_axis)
title('Ramp and 10Hz signal','fontsize',4)
xlabel('Time(sec)','fontsize',3)
ylabel('Level','fontsize',3)
count=4
for i=1:count
Vir_signal(1,1+(i-1)*N*t_range:N*i*t_range)=Signal(1,:)
end
scf(11);plot(0:1/N:count*t_range-1/N,Vir_signal)
title('Virtual analysis signal','fontsize',4)
xlabel('Time(sec)','fontsize',3)
ylabel('Level','fontsize',3)
// HPF を掛ける
s=poly(0,'s')
dt=1/N
dires=[]
cutoffrate=1/N // サンプリング周波数に対する HPF のカットオフ周波数比
fc=1/dt*cutoffrate // カットオフ周波数
HPF=1-fc/(fc+s/2/%pi) // HPF 伝達関数式
HPFsys=syslin('c',HPF);
clf(20);scf(20);clf;bode(HPFsys,1,1e3); // HPF ボーデ線図
title('HPF bode','fontsize',4)
// t=0:dt:dt/cutoffrate
scf(30);clf
HPF_signal=csim(Signal,t,HPFsys) // HPF を通したあとの信号
plot(t,HPF_signal)
title('Signal with HPF','fontsize',4)
xlabel('Time(sec)','fontsize',3)
ylabel('Level','fontsize',3)
// ここまでで波形を発生
fmax=50 // FFT の表示周波数範囲
f0_mix=fft(S_sin(N+1:N*t_range)) // 関数 fft を呼び出す。
M0_mix=abs(f0_mix)*2/N // 振幅を計算
f1_mix=fft(Signal(N+1:N*t_range)) // 関数 fft を呼び出す。
M1_mix=abs(f1_mix)*2/N // 振幅を計算
f2_mix=fft(HPF_signal(N+1:N*(t_range))) // 関数 fft を呼び出す。
M2_mix=abs(f2_mix)*2/N // 振幅を計算
ft_axis=[0,fmax,-1.0,1.0]
// scf(40);clf;plot(0:N-1,M_mix)
// mtlb_axis(ft_axis)
ftdb_axis=[0,fmax,-100,0]
scf(50);clf;plot(0:fmax,20*log10(M0_mix(1:fmax+1)),0:fmax,20*log10(M1_mix(1:fmax+1)),0:fmax,20*log10(M2_mix(1:fmax+1)))
title('Frequency signal level','fontsize',4)
xlabel('Frequency(Hz)','fontsize',3)
ylabel('Level(dB)','fontsize',3)
//xtitle('Frequency sugnal level','Frequency(Hz)','Level(dB)')
legend('only 10Hz','With ramp signal','With ramp and HPF')
mtlb_axis(ftdb_axis)
結果は以下の通りで、手製でやったものとほとんど差がありません。先に書いたように 0Hz にレベルが乗っています。
波形や HPF のボーデ線図は省略しています。
エイリアシングも同様な結果になると思われるので、省略と云うことにします。
Excel:「ところで前から言っていた測定期間内で端数になる周波数はどうすんの？」
画蔵：「それな。新しいイベントボスぐらいの攻撃力があるんで、助けを呼ばなくてはいけない」
FFT 使い手 Level 5 にありがちなこと（ランプ関数の影響を少なくするには）
http://blogs.yahoo.co.jp/susanoo2001_hero/12412381.html
では、手製のスペクトラムアナライザでランプ関数の影響を少なくするための実験を行ってみました。
今度は scilab の FFT で同じことをやってみます。
スクリプトはこれです。
FFT の区間は 1.0 ～ 2.0 秒です。 10Hz 以外は DC が乗ってしまいます。まあこれは結果の 0Hz を無視すると云うことでいいでしょう。
xdel(winsid()) //ウィンドをクリア
clear // Data をクリア
t_range=2 // 測定時間幅
N=2048 // サンプリング周波数
t=0:1/N:t_range-1/N // 時間内のサンプル数
fig_axis=[0,2,-2.0,2.0]
fs=10 // 正弦波の周波数
amp_sin=0.1 // 正弦波の振幅
Ramp(1:N*t_range)=[-(N*t_range)/2+0.5:1:(N*t_range)/2-0.5]/N*2 // ランプ関数
phi=100/N // 正弦波の位相
S_sin(1,1:N*t_range)=sin(2*%pi*(t+phi)*fs)*amp_sin // 本来の正弦波
Signal=Ramp+S_sin
scf(10);clf;plot(t,Signal)
mtlb_axis(fig_axis)
title('Ramp and 10Hz signal','fontsize',4)
xlabel('Time(sec)','fontsize',3)
ylabel('Level','fontsize',3)
count=4
for i=1:count
Vir_signal(1,1+(i-1)*N*t_range:N*i*t_range)=Signal(1,:)
end
scf(11);plot(0:1/N:count*t_range-1/N,Vir_signal)
title('Virtual analysis signal','fontsize',4)
xlabel('Time(sec)','fontsize',3)
ylabel('Level','fontsize',3)
// HPF を掛ける
s=poly(0,'s')
dt=1/N
dires=[]
cutoffrate=1/N // サンプリング周波数に対する HPF のカットオフ周波数比
fc=1/dt*cutoffrate // カットオフ周波数
HPF=1-fc/(fc+s/2/%pi) // HPF 伝達関数式
HPFsys=syslin('c',HPF);
clf(20);scf(20);clf;bode(HPFsys,1,1e3); // HPF ボーデ線図
title('HPF bode','fontsize',4)
// t=0:dt:dt/cutoffrate
scf(30);clf
HPF_signal=csim(Signal,t,HPFsys) // HPF を通したあとの信号
plot(t,HPF_signal)
title('Signal with HPF','fontsize',4)
xlabel('Time(sec)','fontsize',3)
ylabel('Level','fontsize',3)
// ここまでで波形を発生
fmax=50 // FFT の表示周波数範囲
f0_mix=fft(S_sin(N+1:N*t_range)) // 関数 fft を呼び出す。
M0_mix=abs(f0_mix)*2/N // 振幅を計算
f1_mix=fft(Signal(N+1:N*t_range)) // 関数 fft を呼び出す。
M1_mix=abs(f1_mix)*2/N // 振幅を計算
f2_mix=fft(HPF_signal(N+1:N*(t_range))) // 関数 fft を呼び出す。
M2_mix=abs(f2_mix)*2/N // 振幅を計算
ft_axis=[0,fmax,-1.0,1.0]
// scf(40);clf;plot(0:N-1,M_mix)
// mtlb_axis(ft_axis)
ftdb_axis=[0,fmax,-100,0]
scf(50);clf;plot(0:fmax,20*log10(M0_mix(1:fmax+1)),0:fmax,20*log10(M1_mix(1:fmax+1)),0:fmax,20*log10(M2_mix(1:fmax+1)))
title('Frequency signal level','fontsize',4)
xlabel('Frequency(Hz)','fontsize',3)
ylabel('Level(dB)','fontsize',3)
//xtitle('Frequency sugnal level','Frequency(Hz)','Level(dB)')
legend('only 10Hz','With ramp signal','With ramp and HPF')
mtlb_axis(ftdb_axis)
結果は以下の通りで、手製でやったものとほとんど差がありません。先に書いたように 0Hz にレベルが乗っています。
波形や HPF のボーデ線図は省略しています。
エイリアシングも同様な結果になると思われるので、省略と云うことにします。
Excel:「ところで前から言っていた測定期間内で端数になる周波数はどうすんの？」
画蔵：「それな。新しいイベントボスぐらいの攻撃力があるんで、助けを呼ばなくてはいけない」