FFT（Bluestein 和 Cooley-Tukey ）......开始时的意外高峰答案

【问题标题】：FFT (Bluestein and Cooley-Tukey )... unexpected peak at the beginingFFT（Bluestein 和 Cooley-Tukey ）......开始时的意外高峰
【发布时间】：2018-02-06 08:03:30
【问题描述】：

几天前，我发布了我的 FFT 测试工具，但出现了一些问题，所以这是我的“重制”帖子，问题几乎相同，但代码不同。

问题：我正在尝试制作一个工具来读取一些传感器值并执行这些操作。我实现了一个普通的 DFT、一个 FFT（Bluestein）和一个 FFT（Cooley-Tukey）。一切正常，逆变换显示应有的输入。唯一的事情是，当我绘制频率/幅度图时，我在开始时有一个巨大的峰值（不仅像 DC（0Hz）分量应该那样在 0 处！）。

对于 Bluestein 和 DFT，我使用完整输入，对于 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,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 位 ZeroPadding FFT 时

【问题讨论】：

这可能只是“涂抹”的 DC 分量 - 涂抹是由于 spectral leakage 造成的，因为您可能没有应用 window function。
你能显示生成输入数据的代码吗？
我添加了如何生成输入数据
Emm.. 所以Complex 结构包含 x 坐标（斜线）作为实部和一些数据作为虚部？？
sry，但是“斜线”是什么意思？复杂结构在 x 中包含连续区间（newinterval[i] = lastinterval + interval）作为实部，在 Y 中作为虚部包含 y 坐标中正弦波的幅度

标签： c# math graph fft dft

【解决方案1】：

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

【讨论】：

我试过了，它非常接近解决方案。我不得不将真实的部分归零，现在输出是正确的。