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