ST言語の勉強

業務でPLCを用いた製品計測の簡易化・コストダウンに取り組んでおり、その関係でST言語について勉強し始めています。
ST言語はC言語ライクに数値を扱えるので、外部から取り込んだデータを数値計算、数値判定させる用途に適しています。
従来はPCやテスタに取り付けたAD変換ボードでセンサの出力を取り込み、PC上で数値計算・良否判定してPLCに良否判定結果を出力というフローでしたが、これを活用すればすべてをPLC上でクローズさせることができるので構成の簡易化とコストダウンに大いに効果がありそう。
ラズパイをCodesysでPLC化して、ブレッドボード上でいろいろと試行しています。
元々、ラダーはほぼ初心者なのですが、Cはそれなりに経験を積んでいますので、STについてはすっと理解できました。
簡単な自己保持回路等はラダーで組むほうが簡易ですが、少し複雑な処理はSTのほうが見やすく書ける上に再利用が容易なようです。
回路の骨子（メインフロー）はラダーで記述し、個々のブロックは必要に応じてST（インラインST、FB）を組み合わせて記述するのが効率よさそうです。
例　X001のプッシュ毎にY001のON/OFFを切り替える回路（自己保持付き）
ラダー

ST

FFTプログラムの作成（４）

FFTのピーク位置について調べてみました。
FFT実施時にサンプリング個数が2のべき乗に足りない場合は元データを0詰めして2のべき乗にする、というテクニックがありますのでデータ数65536の50001以降のデータを強制的に0にして再度FFT実施しました。
自作、MathNet.Nemricsともにピーク位置は16次に変化ありません。

 続いて強制0ではなく、すべてデータで埋めた50000のサンプリング数でFFTします。2のべき乗ではないので自作プログラムではFFTできませんが、MathNet.NemricsではFFT可能です。
以下がFFT結果です。ピークの位置が12次に移動しました。

どうやらMathNet.Nemrics内部でFFTの手法を切り替えてますね。

調べたところChirp z-変換アルゴリズム（Bluestein FFT）という手法を用いれば2のべき乗でなくともFFTできるようなので、それを使っているのかもしれません。

今回の結果から以下のように理解しました（個人の感想です！）。

情報が不明な信号を解析する場合は0埋めFFTで良いと思いますが、基準の信号波長が既知でその強度を調査したい、等の用途には0埋めは不向きな場合があるようです。
といいますか、65536のデータのうち1/4近くの15536個も0詰めしているのが駄目なんでしょうね。既知の信号を解析する場合は、できるだけ2のべき乗に近しいサンプリング数になるようサンプリングレートを決める必要がある、という当たり前の結論なのか。
いちおう納得したので2のべき乗でないFFTの実装チャレンジしてみます。

FFTプログラムの作成（３）

少し暇ができたので見直し。Complex構造体を使ってみました。ビット反転をもう少し効率化したい。

        public static void FFT2Radix2(ref Complex[] Wave)
        {
            int Number = Wave.Length;
            const double PI = Math.PI;
            int StageCount;
            int BlockCount;
            int NodeCount;
            int StageNumber;
            int BlockNumber;
            int NodeNumber;
            int[] IndexBitChange = new int[Number];

            int TempN1;
            int TempN2;

            Complex Temp1;
            Complex Temp2;
            Complex RotationFactor;
            Complex ButterflyResolution;

            //計算準備
            //累乗数を計算
            StageCount = PrimeFactors(Number).Length;

            //バタフライ演算の繰り返しループ
            for (StageNumber = 0; StageNumber < StageCount; StageNumber++)
            {
                //回転因子の分割数
                BlockCount = (int)(Math.Pow(2, StageNumber));
                //分割した回転因子内の信号数
                NodeCount = Number / BlockCount;
                //回転因子の刻み度数
                ButterflyResolution = new Complex(0, -2 * PI / Number * BlockCount);

                //回転因子の分割数に達するまでループを回す
                for (BlockNumber = 0; BlockNumber < BlockCount; BlockNumber++)
                {
                    //バタフライ演算
                    for (NodeNumber = 0; NodeNumber < NodeCount / 2; NodeNumber++)
                    {
                        //回転因子
                        RotationFactor = Complex.Exp(ButterflyResolution * NodeNumber);
                        //バタフライ演算
                        TempN1 = NodeNumber + NodeCount * BlockNumber;
                        TempN2 = TempN1 + NodeCount / 2;
                        Temp1 = Wave[TempN1];
                        Temp2 = Wave[TempN2];
                        Wave[TempN1] = Temp1 + Temp2;
                        Wave[TempN2] = (Temp1 - Temp2) * RotationFactor;
                    }
                }
            }
            //ビット反転
            for (StageNumber = 0; StageNumber < StageCount; StageNumber++)
            {
                for (int i = 0; i < Math.Pow(2, StageNumber); i++)
                {
                    IndexBitChange[(int)Math.Pow(2, StageNumber) + i] = IndexBitChange[i] +
 　　　　　　　　(int)Math.Pow(2, StageCount - StageNumber - 1);
                }
            }

            for (int i = 0; i < Number; i++)
            {
                if (IndexBitChange[i] > i)
                {
                    Temp1 = Wave[IndexBitChange[i]];
                    Wave[IndexBitChange[i]] = Wave[i];
                    Wave[i] = Temp1;
                }
            }
        }

FFTプログラムの作成（２）

FFTの理解を深めるため、ネットで先達が提示されていたVBのプログラムをC#に移植しました。もっと効率的な書き方できるんだろうけど、理解深めるために愚直に書きます。
実際にN=8=2³の場合で計算順を辿って考えてみるとバタフライ演算の流れが少しずつ理解できました。
最初に大きな回転因子で計算し、それを半分ずつ小さな回転因子に分割していくようなイメージでしょうか？
実際にサンプルデータを作ってみて検算してみました。
        public static void FFT2Radix(double[] InputWave, out double[] OutputWaveReal, out double[] OutputWaveImaginary)
        {
            int Number = InputWave.Length;
            OutputWaveReal = new double[Number];
            OutputWaveImaginary = new double[Number];

            int StageCount;
            int BlockCount;
            int NodeCount;
            int StageNumber;
            int BlockNumber;
            int NodeNumber;
            int[] index = new int[Number];

            int TempN1;
            int TempN2;

            double ButterflyResolution;
            int[] IndexBitChange = new int[Number];
            double Temp1_Real;
            double Temp1_Imaginary;
            double Temp2_Real;
            double Temp2_Imaginary;
            double RotationFactorCos;
            double RotationFactorSin;

            for (int i = 0; i < Number; i++)
            {
                OutputWaveReal[i] = InputWave[i];
                OutputWaveImaginary[i] = 0;
            }

            //計算準備
            //累乗数を計算
            StageCount = PrimeFactors(Number).Length;

            //バタフライ演算の繰り返しループ
            for (StageNumber = 0; StageNumber < StageCount; StageNumber++)
            {
                //回転因子の分割数
                BlockCount = (int)(Math.Pow(2, StageNumber));
                //分割した回転因子内の信号数
                NodeCount = Number / BlockCount;
                //回転因子の刻み度数
                ButterflyResolution = 2 * PI / Number * BlockCount;

                //回転因子の分割数に達するまでループを回す
                for (BlockNumber = 0; BlockNumber < BlockCount; BlockNumber++)
                {
                    //バタフライ演算
                    for (NodeNumber = 0; NodeNumber < NodeCount / 2; NodeNumber++)
                    {
                        RotationFactorCos = Math.Cos(ButterflyResolution * NodeNumber);
                        RotationFactorSin = Math.Sin(ButterflyResolution * NodeNumber);
                        //バタフライ演算
                        TempN1 = NodeNumber + NodeCount * BlockNumber;
                        TempN2 = TempN1 + NodeCount / 2;
                        Temp1_Real = OutputWaveReal[TempN1];
                        Temp1_Imaginary = OutputWaveImaginary[TempN1];
                        Temp2_Real = OutputWaveReal[TempN2];
                        Temp2_Imaginary = OutputWaveImaginary[TempN2];
                        OutputWaveReal[TempN1] = Temp1_Real + Temp2_Real;
                        OutputWaveImaginary[TempN1] = Temp1_Imaginary + Temp2_Imaginary;
                        OutputWaveReal[TempN2] = (Temp1_Real - Temp2_Real) * RotationFactorCos
                            + (Temp1_Imaginary - Temp2_Imaginary) * RotationFactorSin;
                        OutputWaveImaginary[TempN2] = -(Temp1_Real - Temp2_Real) * RotationFactorSin
                            + (Temp1_Imaginary - Temp2_Imaginary) * RotationFactorCos;
                    }
                }
            }
            //インデックス入れ替え
            for (StageNumber = 0; StageNumber < StageCount; StageNumber++)
            {
                for (int i = 0; i < Math.Pow(2, StageNumber); i++)
                {
                    index[(int)Math.Pow(2, StageNumber) + i] = index[i] + (int)Math.Pow(2, StageCount - StageNumber - 1);
                }
            }

            for (int i = 0; i < Number; i++)
            {
                if (index[i] > i)
                {
                    Temp1_Real = OutputWaveReal[index[i]];
                    Temp1_Imaginary = OutputWaveImaginary[index[i]];
                    OutputWaveReal[index[i]] = OutputWaveReal[i];
                    OutputWaveImaginary[index[i]] = OutputWaveImaginary[i];
                    OutputWaveReal[i] = Temp1_Real;
                    OutputWaveImaginary[i] = Temp1_Imaginary;
                }
            }

MathNet.Numericsの計算結果と照らし合わせて検算しましたが全く同じ結果となり、まずはうまくいったようです。元の波形は以下。
Vector[i] = 3 + 0.5 * Math.Cos(2 * Math.PI * 12.0 * 2.0* (i / 100000))
振幅0.5Vの正弦波にDCが3V乗ったイメージです。
直流成分をサンプル数で割るとほぼ3になるので問題なし。ただし本来なら12次の位置にピークが出ると思っていたら16次の位置に出ています。
理由は何となくわかりましたが、ここまでで疲れたのでまた次回にします。

FFTプログラムの作成（１）

信号波形のスペクトル解析を実施する機会が結構な頻度あり、.NET Framework向けのオープンソースな数値計算用のライブラリを使ってC#で解析プログラム作成してました。

Math.Numerics
https://numerics.mathdotnet.com/

このライブラリは非常に使い勝手いいのですが、今後もオープンソースであり続けるかという点で個人的に不安感あるのと、FFTはデータ数が2の累乗でなければならない（というか、最も計算効率が良い）ということですが、上記ライブラリ使用時は気にせず不定長のデータを入力してます。オシロスコープ上で実施するスペクトル解析の信号強度とほぼ同じ値（誤差±0.2mV以下）が出てくるので実質上は問題ないのですが、そこがブラックボックスなのでどうにも気持ちが悪いのです。

上記の理由より、オリジナルでライブラリ作成を検討しています。

FFTについてアルゴリズムの概要はだいたい理解できたのですが、バタフライ演算をどうプログラムに実装するか、が中々理解できません。
ネットで調べるといろいろな方がバタフライ演算のサンプルプログラムを提示してくれていますが、どうにも腹落ちができてません。
まずは形から入り覚える、ということで先達のプログラムをC#に移植しています。
とりあえず自分の頭の整理と覚えのため記載しました。