FFT (Bluestein和Cooley-Tukey )...始料未及的高峰

Question

几天前，我发布了我的FFT-Test工具，但有一些问题，所以这里是我的“重做”帖子，几乎相同的问题，但不同的代码。

问题:我正在尝试制作一个工具来读取一些传感器的值，并对这些值进行fft。我实现了一个普通的DFT，一个FFT(Bluestein)和一个FFT(Cooley-Tukey)。这一切都很好，并且逆变换显示了它应该显示的输入。唯一的问题是，当我绘制Frequenz/振幅图时，我在一开始就有一个巨大的峰值(而不是像DC (0Hz)分量那样只在0处！)。

对于Bluestein和DFT，我使用complet输入，对于Cooley-Tukey，我使用ZeroPadding计算，或者只是将输入削减到2^n的长度。

例如，这是我的输入，只是一个简单的正弦波。这里输出的是DFT和Bluetstein(看起来是一样的)、FFT(ZeroPadding)和一个DFT zoomed in at 0Hz。这是用八度音阶输入的ouput，这是我期望得到的结果。

这里是我的CT-FFT之前的代码

            public static void TransformRadix2(Complex[] vector, bool inverse)
        {
            // Length variables
            int n = vector.Length;
            int levels = 0;  // compute levels = floor(log2(n))
            for (int temp = n; temp > 1; temp >>= 1)
                levels++;
            if (1 << levels != n)
                throw new ArgumentException("Length is not a power of 2");

            // Trigonometric table
            Complex[] expTable = new Complex[n / 2];
            double coef = 2 * Math.PI / n * (inverse ? 1 : -1);
            for (int i = 0; i < n / 2; i++)
                expTable[i] = Complex.Exp(new Complex(0, i * coef));

            // Bit-reversed addressing permutation
            for (int i = 0; i < n; i++)
            {
                int j = (int)((uint)ReverseBits(i) >> (32 - levels));
                if (j > i)
                {
                    Complex temp = vector[i];
                    vector[i] = vector[j];
                    vector[j] = temp;
                }
            }

            // Cooley-Tukey decimation-in-time radix-2 FFT
            for (int size = 2; size <= n; size *= 2)
            {
                int halfsize = size / 2;
                int tablestep = n / size;
                for (int i = 0; i < n; i += size)
                {
                    for (int j = i, k = 0; j < i + halfsize; j++, k += tablestep)
                    {
                        Complex temp = vector[j + halfsize] * expTable[k];
                        vector[j + halfsize] = vector[j] - temp;
                        vector[j] += temp;
                    }
                }
                if (size == n)  // Prevent overflow in 'size *= 2'
                    break;
            }
        }

这是布鲁斯坦。

            public static void TransformBluestein(Complex[] vector, bool inverse)
        {
            // Find a power-of-2 convolution length m such that m >= n * 2 + 1
            int n = vector.Length;
            if (n >= 0x20000000)
                throw new ArgumentException("Array too large");
            int m = 1;
            while (m < n * 2 + 1)
                m *= 2;

            // Trignometric table
            Complex[] expTable = new Complex[n];
            double coef = Math.PI / n * (inverse ? 1 : -1);
            for (int i = 0; i < n; i++)
            {
                int j = (int)((long)i * i % (n * 2));  // This is more accurate than j = i * i
                expTable[i] = Complex.Exp(new Complex(0, j * coef));
            }

            // Temporary vectors and preprocessing
            Complex[] avector = new Complex[m];
            for (int i = 0; i < n; i++)
                avector[i] = vector[i] * expTable[i];
            Complex[] bvector = new Complex[m];
            bvector[0] = expTable[0];
            for (int i = 1; i < n; i++)
                bvector[i] = bvector[m - i] = Complex.Conjugate(expTable[i]);

            // Convolution
            Complex[] cvector = new Complex[m];
            Convolve(avector, bvector, cvector);

            // Postprocessing
            for (int i = 0; i < n; i++)
                vector[i] = cvector[i] * expTable[i];
        }


        /* 
         * Computes the circular convolution of the given complex vectors. Each vector's length must be the same.
         */
        public static void Convolve(Complex[] xvector, Complex[] yvector, Complex[] outvector)
        {
            int n = xvector.Length;
            if (n != yvector.Length || n != outvector.Length)
                throw new ArgumentException("Mismatched lengths");
            xvector = (Complex[])xvector.Clone();
            yvector = (Complex[])yvector.Clone();
            Transform(xvector, false);
            Transform(yvector, false);
            for (int i = 0; i < n; i++)
                xvector[i] *= yvector[i];
            Transform(xvector, true);
            for (int i = 0; i < n; i++)  // Scaling (because this FFT implementation omits it)
                outvector[i] = xvector[i] / n;
        }

        private static int ReverseBits(int val)
        {
            int result = 0;
            for (int i = 0; i < 32; i++, val >>= 1)
                result = (result << 1) | (val & 1);
            return result;
        }

这是每个周期的输入，间隔为0,00001953125。

0
 0.03141
 0.06279
 0.09411
 0.12533
 0.15643
0.18738
 0.21814
 0.24869
 0.27899
 0.30902
 0.33874
 0.36812
 0.39715
 0.42578
 0.45399
 0.48175
 0.50904
 0.53583
0.56208
 0.58779
 0.61291
 0.63742
 0.66131
 0.68455
0.70711
 0.72897
 0.75011
 0.77051
 0.79016
 0.80902
 0.82708
 0.84433
 0.86074
 0.87631
 0.89101
 0.90483
 0.91775
 0.92978
 0.94088
 0.95106
 0.96029
 0.96858
 0.97592
0.98229
 0.98769
 0.99211
 0.99556
 0.99803
 0.99951
1
 0.99951
   0.99803
   0.99556
   0.99211
   0.98769
   0.98229
   0.97592
   0.96858
   0.96029
   0.95106
   0.94088
   0.92978
  0.91775
   0.90483
   0.89101
   0.87631
   0.86074
   0.84433
  0.82708
  0.80902
  0.79016
   0.77051
   0.75011
   0.72897
   0.70711
   0.68455
   0.66131
   0.63742
   0.61291
   0.58779
   0.56208
   0.53583
   0.50904
   0.48175
   0.45399
   0.42578
   0.39715
  0.36812
   0.33874
   0.30902
   0.27899
   0.24869
   0.21814
  0.18738
   0.15643
   0.12533
   0.09411
   0.06279
   0.03141
   0
   -0.03141
   -0.06279
   -0.09411
   -0.12533
   -0.15643
   -0.18738
  -0.21814
   -0.24869
   -0.27899
   -0.30902
   -0.33874
   -0.36812
-0.39715
   -0.42578
   -0.45399
   -0.48175
   -0.50904
   -0.53583
  -0.56208
   -0.58779
   -0.61291
   -0.63742
   -0.66131
   -0.68455
   -0.70711
   -0.72897
   -0.75011
   -0.77051
   -0.79016
   -0.80902
   -0.82708
  -0.84433
   -0.86074
   -0.87631
   -0.89101
   -0.90483
   -0.91775
  -0.92978
   -0.94088
   -0.95106
   -0.96029
   -0.96858
   -0.97592
   -0.98229
   -0.98769
   -0.99211
   -0.99556
   -0.99803
   -0.99951
   -1
  -0.99951
   -0.99803
   -0.99556
   -0.99211
   -0.98769
   -0.98229
  -0.97592
   -0.96858
   -0.96029
   -0.95106
   -0.94088
   -0.92978
  -0.91775
   -0.90483
   -0.89101
   -0.87631
   -0.86074
   -0.84433
   -0.82708
   -0.80902
   -0.79016
   -0.77051
   -0.75011
   -0.72897
   -0.70711
  -0.68455
   -0.66131
   -0.63742
   -0.61291
   -0.58779
   -0.56208
  -0.53583
   -0.50904
   -0.48175
   -0.45399
   -0.42578
   -0.39715
   -0.36812
   -0.33874
   -0.30902
   -0.27899
   -0.24869
   -0.21814
   -0.18738
   -0.15643
   -0.12533
   -0.09411
   -0.06279
   -0.03141
   0

我不知道这个高峰会是什么样子。我唯一能说的是，当我的区间越大，开头的峰值也越大: interval 0,00001953125；0,0001953125；0,001953125

下面是我如何创建fft的输入。只有我显示的输入和间隔....

                    float UsedInterval = float.Parse(Interval.Replace(",", "."), CultureInfo.InvariantCulture);
                float x = 0; //Intervall starting at 0
                foreach (var line in input)
                {
                    if (string.IsNullOrWhiteSpace(line)) continue; //delete empty lines

                    x += UsedInterval;  //add intervall
                    double y = double.Parse(line.Replace(",", "."), CultureInfo.InvariantCulture);

                    arr.Add(new Vector2XY(x, y));
                }

我希望你能帮助我解决我的问题，谢谢你的关注。

当我用一个周期上的高位数据和间隔0,00001953125 this is my ouput进行256位FFT时

Answer 1

你应该用Sin(coeff*i)填充复数输入的实部，用零-纯实输入填充虚部，这样你就会得到正确的输出，在零附近没有大的峰值-这是由于当前实部的线性分量。