FFTを使った、直線たたき込み
下の low pass filter を、+5KHzだけ、複素周波数遷移した
1601tap low pass filter ( linear convolution using FFT )
1.FFTを使った 直線たたき込み ...「.直線たたみ込み」 でしょっ! 又、間違えた (^_^;;
filterを作るのに、窓関数法を使いました。
正規化された遷移周波数 delta_frequencyから、逆算しました。
窓関数法によるfilter の設計は、
「ディジタル サウンド処理入門 青木直史先生 CQ出版社」
が、超、詳しくて、解りやすいです。
以下は、私のプログラムの一部です。
namespace
myDirectSoundCapture1{
public partial class MainForm : Form
{
private double fc = 0; //cut off frequency
private int fe = 0; //edge frequency
private int filter_size;
private double delta_frequency_current;
private double[] taps_current ;
//窓関数法による、fir filterの作成
public double[] get_LPFir_coeff( int size, int f )
{
filter_size = size;
fe = f; //edge frequency
// use window
switch (WindowFunc_comboBox.SelectedIndex)
{
case 0: delta_frequency_current = (double)(0.9 / filter_size); break; //rectangular
case 1: delta_frequency_current = (double)(3.1 / filter_size); break; //Hanning
case 2: delta_frequency_current = (double)(3.3 / filter_size); break; //Hamming
case 3: delta_frequency_current = (double)(5.5 / filter_size); break; //Blackman
default: delta_frequency_current = (double)(3.3 / filter_size); break; //default
}
delta_frequency_current *= my_waveformat.SamplesPerSecond;//de_normalize
fc = fe + delta_frequency_current / 2.0; // not normalized
//normalize
fc /= my_waveformat.SamplesPerSecond;
taps_current = new double[filter_size*4];
int K = 0;
for (int i = -filter_size*2; i < filter_size*2; i++) //filterを、ずらさないと、インパルス応答が、半分しか出ない。
{
if (i == 0)
taps_current[K] = 2.0 * fc;
else
{
taps_current[K] = 2.0 * fc * Math.Sin(2 * Math.PI * fc * i) / (2 * Math.PI * fc * i);
}
K++;
}
do_window( taps_current.Length, taps_current, taps_current );
return taps_current;
}
.etc..........................
作った関数を、使う所です、一部抜粋。
cos_FFT(.double[] coef..)は、実係数( coef )を、複素周波数遷移させる関数。
2.使用する窓関数による差
Rectangular窓での出力を0dBとします。
+5KHz comlex shift
下図は、拡大したものです。
edge freq.より、+46.5Hzで、約-80dB、落ちるが、出力も、rectangular窓に比べて、
1/8 位 (-18dB)、落ちる。(-6dB=1/2、-12dB=1/4)
Hamming 窓
Rectangularと比べて、-10dB down , 出力が、1/4近くに、下がる。( -6dB = 1/2 )
+5KHz shift ( using complex )
Blackman 窓
+5KHz shift
3.気になる所
9.5KHz辺りの出力を計測している時、filterの通過帯域にも、大きな出力が見られる事。
-46dB位下では、ありまするが.....
通過帯域の出力に比べて、1/200位(40dB=1/100、6dB=その半分)
9.5KHzの出力自体は、-80dBされるのでは、ありますが、その影響は、通過帯域にも出現し、
結局 -46dBcが通過帯域に、出力されています。
ひどい時には、以下のように、なります。
1.5KHz成分は、-90dBされるが、その影響は、通過帯域で、-38dBc位の雑音?となって、観測される。
由々しき問題を、発見いたしましたわ.....涙、涙.....
もしかしたら、通常の 時間領域の、convolutionでは、出現しないかも.....と、期待して
時間領域での、通常の convolutionを行ってみました。
あかんがな....やっぱり、取れませんわ......
時間領域でのconvolutionは、原理そのままで、行いました。
以下は、そのプログラムの一部です。
尚、処理に要する時間は、FFTを使った直線たたみ込みが、圧倒的に、有利でした。
4.偶数tapの方が、結果がよい
2049tap Hanning窓使用の、窓関数法によるLPF + freq.shift( +5KHz )
通過帯域が、いびつです。
2048tap Hanning LPF + freq.shift( +5KHz )
2048tap の偶数の方が、通過帯域が平坦で、素直です。
窓関数法によるfilterでは、[フィルタのサイズは、奇数にする必要がある」らしいのですが
偶数の方が、結果が、よろしい..... 不思議やな.....
1000 tapでも、かなりよい。遷移周波数も、75Hz位と、いいですね。
時間的な、delayも、大分ましになると、思います。1000 tapでも、充分かもしれませんね。
先ほどの、通過帯域の 雑音は、
入力レベルを調節しても、firの係数のscalingを調整しても、変わりませんでした。
通過帯域のトップレベルと、雑音のレベルが、共に、平行移動するだけでした。
軽減する方法は、ないか?
通過帯域外の信号が、弱いなら、問題ありません。その信号の強さに応じて、通過帯域内の雑音も
減少します。
帯域外の信号は、強力に減衰してくれるのに、帯域内への、影響が大きいとは.....知らなかった。
そう言う事やったんか....実験して、初めて気づきましたわ。
サンプリングが、16bitではなくて、24bitでは、どうなんやろ....段階がきめ細かくなるだけ、やろなあ....
24bitを試すには、プログラムが、めちゃ難しくなる...C#が使えない...今ん所.......
今回は、これ(16bit)で、実験を続けてみましょうか....
H.18.10.19