今回は FFT だけに関係するわけではないですが、アナログ信号をデジタル信号に変換する際に必ずついて回るエイリアシングについて、再度確認してみようと思います。
一度はここで触れています。
サンプリングによるスペクトラムの変化（エイリアシングの確認）
http://blogs.yahoo.co.jp/susanoo2001_hero/8968481.html
さて、このエイリアシングの影響を除くためにサンプリングというか普通は AD コンバータにアナログ信号を入力する前にフィルタを入れるわけですが、ナイキスト周波数（サンプリング周波数の半分）の直前で急峻な減衰特性を持つ LPF の実現が大変です。
今はそんなことはないと思いますが、昔のデジタルオシロでは上述の前置 LPF（アンチエイリアシングフィルタ）が入っているにはいるのですが、タイムスケールを変えた時、AD コンバータのサンプリング周波数だけを変化させて、前置 LPF の特性を変化させていなかっったようで、高周波ノイズの多い信号を低い周波数のみ観測しようとすると、奇妙な波形が乗り込んできて悩まされた記憶があります。
多分、今はそんなことないでしょう。とりあえず超高速のサンプリング周波数で AD 変換してながら必要なタイムスケールに間引く前に間引く周波数に合わせて強力な FIR フィルタで減衰させてしまえば良いので、そういった問題は起きないわけです。
以上の内容はここでもちょっと触れています。
オーバーサンプリングをしてみる＞エイリアシングを防ぐには
http://blogs.yahoo.co.jp/susanoo2001_hero/8974066.html?type=folderlist
ナイキスト周波数の直前まで本当に再現出来るのか？＞標本化定理
http://blogs.yahoo.co.jp/susanoo2001_hero/8998742.html?type=folderlist
FIR デジタルフィルタで理想 LPF もどきを作ってみる
http://blogs.yahoo.co.jp/susanoo2001_hero/9168591.html?type=folderlist
なんだ、過去の書き込みの紹介だらけではないか、と思われてしまいそうなのですが、実際に使うにあたっていつもいつも高速の AD コンバータがあるわけでもないし、急峻な FIR を実現できるほどのロジック回路規模があるわけでもないし、というのが実態だしまた測定精度は組み込み機器ならこの程度、といったこともあると思うので、高い周波数に対する落としどころを探ってみようというわけです。
今回のスクリプトのイメージは、前回までのスクリプトで使われていたサンプリング周波数を 10倍にして、これをアナログ信号と見なして、ここに今までと同じサンプリング周波数でのナイキスト周波数を越えるノイズが混入しているとし、シンプルな一次 LPF でノイズ除去をしてから見たい周波数にどのような影響を与えるか、ということになります。
xdel(winsid()) //ウィンドをクリア
clear // Data をクリア
over_range=10 // オーバサンプリング比
N=2000 // サンプリング周波数: 実際のサンプリング周波数＝over_range x N
No=N*over_range
to=0:1/No:1-1/No // 時間内のサンプル数
t=0:1/N:1-1/N //
fig_axis=[0,2,-2.0,2.0]
fs=100 // 正弦波の周波数
amp_s=1 // 正弦波の振幅 0dB
fn=1890 // ノイズの周波数
amp_n=0.5 // ノイズの振幅 -6dB
phi_n=100/N // ノイズの位相
S_sin(1,1:No)=sin(2*%pi*to*fs)*amp_s // 本来の正弦波
N_sin(1,1:No)=sin(2*%pi*(to+phi_n)*fn)*amp_n // ノイズの正弦波
Mix_signal1=S_sin+N_sin // 合成された信号
f_min=1
f_max=N
for f1=f_min:f_max
x=sin(2*%pi*to*(f1)) // 正弦波
y=cos(2*%pi*to*(f1)) // 余弦波
combo_sin1=x.*Mix_signal1
combo_cos1=y.*Mix_signal1
sin_level1(1,f1)=sum(combo_sin1)./No*2
cos_level1(1,f1)=sum(combo_cos1)./No*2
level1(1,f1)=((sin_level1(f1))^2+(cos_level1(f1))^2)^0.5
end
ftdb_axis=[f_min,f_max,-100,0]
scf(20);clf;plot(f_min:f_max,20*log10(level1(f_min:f_max)))
mtlb_axis(ftdb_axis)
title('Original signal spectrum','fontsize',4)
xlabel('(Frequnecy(Hz))','fontsize',3)
ylabel('Level(dB)','fontsize',3)
Mix_signal=Mix_signal1(1:10:No) // ダウンサンプリング
t_axis=[0,0.2,-2.0,2.0]
scf(10);
subplot(2,1,2)
plot(t,Mix_signal) // 入力された信号波形
title('Down sampling signal','fontsize',4)
mtlb_axis(t_axis)
xlabel('Time(sec)','fontsize',3)
ylabel('Level','fontsize',3)
subplot(2,1,1)
plot(to,Mix_signal1) // 単純にダウンサンプリングした時の波形
title('Original signal','fontsize',4)
mtlb_axis(t_axis)
xlabel('Time(sec)','fontsize',3)
ylabel('Level','fontsize',3)
f_min=1
f_max=2000
f_count=f_max-f_min+1
for f1=f_min:f_max
x=sin(2*%pi*t*(f1)) // 正弦波
y=cos(2*%pi*t*(f1)) // 余弦波
combo_sin=x.*Mix_signal
combo_cos=y.*Mix_signal
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
ftdb_axis=[f_min,f_max,-100,0]
scf(30);clf;plot(f_min:f_max,20*log10(level(f_min:f_max))) // 単純にダウンサンプリングした時のスペクトラム
mtlb_axis(ftdb_axis)
title('Down sampled signal spectrum','fontsize',4)
xlabel('(Frequnecy(Hz))','fontsize',3)
ylabel('Level(dB)','fontsize',3)
s=poly(0,'s')
dt=1/No
dires=[]
cutoffrate=20/N // サンプリング周波数に対する HPF のカットオフ周波数比。ここでは 200Hz
fc=1/dt*cutoffrate // カットオフ周波数
LPF=fc/(fc+s/2/%pi) // LPF 伝達関数式
LPFsys=syslin('c',LPF);
scf(40);clf;bode(LPFsys,10,1e4); // HPF ボーデ線図
title('LPF bode','fontsize',4)
LPF_signal1=csim(Mix_signal1,to,LPFsys) // LPF を通したあとの信号
scf(50)
subplot(2,1,1)
plot(to(1:4000),LPF_signal1(1:4000))
title('Signal with LPF','fontsize',4)
xlabel('Time(sec)','fontsize',3)
ylabel('Level','fontsize',3)
mtlb_axis(t_axis)
LPF_signal=LPF_signal1(1:10:No) // LPF 後のダウンサンプリング
subplot(2,1,2)
plot(t,LPF_signal)
title('Down ampling signal with LPF','fontsize',4)
xlabel('Time(sec)','fontsize',3)
ylabel('Level','fontsize',3)
mtlb_axis(t_axis)
f_min=1
f_max=200
for f1=f_min:f_max
x=sin(2*%pi*t*(f1)) // 正弦波
y=cos(2*%pi*t*(f1)) // 余弦波
combo_sin2=x.*LPF_signal
combo_cos2=y.*LPF_signal
sin_level2(1,f1)=sum(combo_sin2)./N*2
cos_level2(1,f1)=sum(combo_cos2)./N*2
level2(1,f1)=((sin_level2(f1))^2+(cos_level2(f1))^2)^0.5
end
ftdb_axis=[f_min,f_max,-200,0]
scf(60);
plot(f_min:f_max,20*log10(level(f_min:f_max)),f_min:f_max,20*log10(level2(f_min:f_max)))
title('Spectrum','fontsize',4)
xlabel('frequency(Hz)','fontsize',3)
ylabel('Level(dB)','fontsize',3)
mtlb_axis(ftdb_axis)
legend('Without LPF','With LPF')
disp(20*log10(level(100))) // 観測したい信号レベル
disp(20*log10(level(110))) // エイリアシングによる信号レベル
disp(20*log10(level2(100))) // LPF を通してから観測した信号レベル
disp(20*log10(level2(110))) // LPF を通してから観測したエイリアシングによる信号レベル
かなり長くなってしまいました。それと実際にスクリプトを実行すると、DFT のところで周波数範囲が広いのでちょっと時間が掛かります。
ということで、結果のグラフなどの解説は次回にしたいと思います。
最後に表示される値のだけのせておきます。
観測したい 100Hz のレベル： - 9.643D-16(dB) （ほとんど 0dB）
1890Hz がエイリアシングによって 110Hz に現れたレベル： - 6.0205999(dB)
LPF を通して後の 100Hz のレベル： - 0.9707148(dB)
LPF を通した後の 110Hz のレベル： - 25.736611(dB)
一度はここで触れています。
サンプリングによるスペクトラムの変化（エイリアシングの確認）
http://blogs.yahoo.co.jp/susanoo2001_hero/8968481.html
さて、このエイリアシングの影響を除くためにサンプリングというか普通は AD コンバータにアナログ信号を入力する前にフィルタを入れるわけですが、ナイキスト周波数（サンプリング周波数の半分）の直前で急峻な減衰特性を持つ LPF の実現が大変です。
今はそんなことはないと思いますが、昔のデジタルオシロでは上述の前置 LPF（アンチエイリアシングフィルタ）が入っているにはいるのですが、タイムスケールを変えた時、AD コンバータのサンプリング周波数だけを変化させて、前置 LPF の特性を変化させていなかっったようで、高周波ノイズの多い信号を低い周波数のみ観測しようとすると、奇妙な波形が乗り込んできて悩まされた記憶があります。
多分、今はそんなことないでしょう。とりあえず超高速のサンプリング周波数で AD 変換してながら必要なタイムスケールに間引く前に間引く周波数に合わせて強力な FIR フィルタで減衰させてしまえば良いので、そういった問題は起きないわけです。
以上の内容はここでもちょっと触れています。
オーバーサンプリングをしてみる＞エイリアシングを防ぐには
http://blogs.yahoo.co.jp/susanoo2001_hero/8974066.html?type=folderlist
ナイキスト周波数の直前まで本当に再現出来るのか？＞標本化定理
http://blogs.yahoo.co.jp/susanoo2001_hero/8998742.html?type=folderlist
FIR デジタルフィルタで理想 LPF もどきを作ってみる
http://blogs.yahoo.co.jp/susanoo2001_hero/9168591.html?type=folderlist
なんだ、過去の書き込みの紹介だらけではないか、と思われてしまいそうなのですが、実際に使うにあたっていつもいつも高速の AD コンバータがあるわけでもないし、急峻な FIR を実現できるほどのロジック回路規模があるわけでもないし、というのが実態だしまた測定精度は組み込み機器ならこの程度、といったこともあると思うので、高い周波数に対する落としどころを探ってみようというわけです。
今回のスクリプトのイメージは、前回までのスクリプトで使われていたサンプリング周波数を 10倍にして、これをアナログ信号と見なして、ここに今までと同じサンプリング周波数でのナイキスト周波数を越えるノイズが混入しているとし、シンプルな一次 LPF でノイズ除去をしてから見たい周波数にどのような影響を与えるか、ということになります。
xdel(winsid()) //ウィンドをクリア
clear // Data をクリア
over_range=10 // オーバサンプリング比
N=2000 // サンプリング周波数: 実際のサンプリング周波数＝over_range x N
No=N*over_range
to=0:1/No:1-1/No // 時間内のサンプル数
t=0:1/N:1-1/N //
fig_axis=[0,2,-2.0,2.0]
fs=100 // 正弦波の周波数
amp_s=1 // 正弦波の振幅 0dB
fn=1890 // ノイズの周波数
amp_n=0.5 // ノイズの振幅 -6dB
phi_n=100/N // ノイズの位相
S_sin(1,1:No)=sin(2*%pi*to*fs)*amp_s // 本来の正弦波
N_sin(1,1:No)=sin(2*%pi*(to+phi_n)*fn)*amp_n // ノイズの正弦波
Mix_signal1=S_sin+N_sin // 合成された信号
f_min=1
f_max=N
for f1=f_min:f_max
x=sin(2*%pi*to*(f1)) // 正弦波
y=cos(2*%pi*to*(f1)) // 余弦波
combo_sin1=x.*Mix_signal1
combo_cos1=y.*Mix_signal1
sin_level1(1,f1)=sum(combo_sin1)./No*2
cos_level1(1,f1)=sum(combo_cos1)./No*2
level1(1,f1)=((sin_level1(f1))^2+(cos_level1(f1))^2)^0.5
end
ftdb_axis=[f_min,f_max,-100,0]
scf(20);clf;plot(f_min:f_max,20*log10(level1(f_min:f_max)))
mtlb_axis(ftdb_axis)
title('Original signal spectrum','fontsize',4)
xlabel('(Frequnecy(Hz))','fontsize',3)
ylabel('Level(dB)','fontsize',3)
Mix_signal=Mix_signal1(1:10:No) // ダウンサンプリング
t_axis=[0,0.2,-2.0,2.0]
scf(10);
subplot(2,1,2)
plot(t,Mix_signal) // 入力された信号波形
title('Down sampling signal','fontsize',4)
mtlb_axis(t_axis)
xlabel('Time(sec)','fontsize',3)
ylabel('Level','fontsize',3)
subplot(2,1,1)
plot(to,Mix_signal1) // 単純にダウンサンプリングした時の波形
title('Original signal','fontsize',4)
mtlb_axis(t_axis)
xlabel('Time(sec)','fontsize',3)
ylabel('Level','fontsize',3)
f_min=1
f_max=2000
f_count=f_max-f_min+1
for f1=f_min:f_max
x=sin(2*%pi*t*(f1)) // 正弦波
y=cos(2*%pi*t*(f1)) // 余弦波
combo_sin=x.*Mix_signal
combo_cos=y.*Mix_signal
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
ftdb_axis=[f_min,f_max,-100,0]
scf(30);clf;plot(f_min:f_max,20*log10(level(f_min:f_max))) // 単純にダウンサンプリングした時のスペクトラム
mtlb_axis(ftdb_axis)
title('Down sampled signal spectrum','fontsize',4)
xlabel('(Frequnecy(Hz))','fontsize',3)
ylabel('Level(dB)','fontsize',3)
s=poly(0,'s')
dt=1/No
dires=[]
cutoffrate=20/N // サンプリング周波数に対する HPF のカットオフ周波数比。ここでは 200Hz
fc=1/dt*cutoffrate // カットオフ周波数
LPF=fc/(fc+s/2/%pi) // LPF 伝達関数式
LPFsys=syslin('c',LPF);
scf(40);clf;bode(LPFsys,10,1e4); // HPF ボーデ線図
title('LPF bode','fontsize',4)
LPF_signal1=csim(Mix_signal1,to,LPFsys) // LPF を通したあとの信号
scf(50)
subplot(2,1,1)
plot(to(1:4000),LPF_signal1(1:4000))
title('Signal with LPF','fontsize',4)
xlabel('Time(sec)','fontsize',3)
ylabel('Level','fontsize',3)
mtlb_axis(t_axis)
LPF_signal=LPF_signal1(1:10:No) // LPF 後のダウンサンプリング
subplot(2,1,2)
plot(t,LPF_signal)
title('Down ampling signal with LPF','fontsize',4)
xlabel('Time(sec)','fontsize',3)
ylabel('Level','fontsize',3)
mtlb_axis(t_axis)
f_min=1
f_max=200
for f1=f_min:f_max
x=sin(2*%pi*t*(f1)) // 正弦波
y=cos(2*%pi*t*(f1)) // 余弦波
combo_sin2=x.*LPF_signal
combo_cos2=y.*LPF_signal
sin_level2(1,f1)=sum(combo_sin2)./N*2
cos_level2(1,f1)=sum(combo_cos2)./N*2
level2(1,f1)=((sin_level2(f1))^2+(cos_level2(f1))^2)^0.5
end
ftdb_axis=[f_min,f_max,-200,0]
scf(60);
plot(f_min:f_max,20*log10(level(f_min:f_max)),f_min:f_max,20*log10(level2(f_min:f_max)))
title('Spectrum','fontsize',4)
xlabel('frequency(Hz)','fontsize',3)
ylabel('Level(dB)','fontsize',3)
mtlb_axis(ftdb_axis)
legend('Without LPF','With LPF')
disp(20*log10(level(100))) // 観測したい信号レベル
disp(20*log10(level(110))) // エイリアシングによる信号レベル
disp(20*log10(level2(100))) // LPF を通してから観測した信号レベル
disp(20*log10(level2(110))) // LPF を通してから観測したエイリアシングによる信号レベル
かなり長くなってしまいました。それと実際にスクリプトを実行すると、DFT のところで周波数範囲が広いのでちょっと時間が掛かります。
ということで、結果のグラフなどの解説は次回にしたいと思います。
最後に表示される値のだけのせておきます。
観測したい 100Hz のレベル： - 9.643D-16(dB) （ほとんど 0dB）
1890Hz がエイリアシングによって 110Hz に現れたレベル： - 6.0205999(dB)
LPF を通して後の 100Hz のレベル： - 0.9707148(dB)
LPF を通した後の 110Hz のレベル： - 25.736611(dB)