FFT(Bluestein和Cooley-Tukey)...在开始时意外的高峰 [英] FFT (Bluestein and Cooley-Tukey )... unexpected peak at the begining

查看:244
本文介绍了FFT(Bluestein和Cooley-Tukey)...在开始时意外的高峰的处理方法,对大家解决问题具有一定的参考价值,需要的朋友们下面随着小编来一起学习吧!

问题描述

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



问题:
我正尝试制作一个工具来读取一些传感器值并做这些。我实现了一个正常的DFT,一个FFT(Bluestein)和一个FFT(Cooley-Tukey)。它一切正常,逆转换显示输入应该像它应该。
唯一的是当我绘制Frequenz /幅度图时,我在开始时有一个巨大的峰值(不仅在0像DC(0Hz)分量一样)!



对于Bluestein和DFT,我使用Complet输入,对于Cooley-Tukey,我使用ZeroPadding进行计算,或者将输入剪切为2 ^ n的长度。 b
$ b

例如,这里是我的输入,只是一个简单的正弦波。
这里DFT和Bluetstein的输出(看起来相同)和 FFT(ZeroPadding)和一个 DFT以0Hz放大
这里是八度音的相同输入的输出,我的预计会得到。



在这里,我的代码前CT-FFT

  public static void TransformRadix2(Complex [] vector,bool inverse)
{
//长度变量
int n = vector.Length;
int levels = 0; (int temp = n; temp> 1; temp>> = 1)
levels ++; //计算级别= floor(log2(n))
;
if(1 <抛出新的ArgumentException(长度不是2的幂);

//三角函数表
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));

//位反转寻址置换
for(int i = 0; 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; (int j = i,k = 0; j< i + halfsize; j ++, k + = tablestep)
{
复数temp = vector [j + halfsize] * expTable [k];
vector [j + halfsize] = vector [j] - temp;
vector [j] + = temp;


if(size == n)//防止'size * = 2'溢出
break;
}
}

此处Bluestein。

  public static void TransformBluestein(Complex [] vector,bool inverse)
{
//找到2的幂卷积长度m使得m> = n * 2 + 1
int n = vector.Length;
if(n> = 0x20000000)
抛出新的ArgumentException(Array too large);
int m = 1; $ m $ b而(m m * = 2;

//三角函数表
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)); //这比j = i * i
expTable [i] = Complex.Exp(new Complex(0,j * coef))更精确;
}

//临时向量和预处理
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]);

//卷积
复数[] cvector = new Complex [m];
卷积(avector,bvector,cvector);

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


/ *
*计算给定复矢量的循环卷积。每个矢量的长度必须相同。
* /
public static void Convolve(Complex [] xvector,Complex [] yvector,Complex [] outvector)
{
int n = xvector.Length;
if(n!= yvector.Length || n!= outvector.Length)
抛出新的ArgumentException(不匹配的长度);
xvector =(Complex [])xvector.Clone();
yvector =(Complex [])yvector.Clone();
Transform(xvector,false);
变换(yvector,false);
for(int i = 0; i< n; i ++)
xvector [i] * = yvector [i];
Transform(xvector,true);
for(int i = 0; i outvector [i] = xvector [i] / n;
}

private static int ReverseBits(int val)
{
int result = 0; (int i = 0; i< 32; i ++,val>> = 1)
结果=(结果<< 1) (val& 1);
返回结果;
}

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

  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

我不知道这个高峰可能是什么。我唯一能说的是,当我的间隔更大时,开始的峰值也更大:
间隔 0,00001953125 ; 0,0001953125 ; 0,001953125



这里是我如何创建fft的输入。只是我显示的输入和intervall ....

  float UsedInterval = float.Parse(Interval.Replace(, 。),CultureInfo.InvariantCulture); 
float x = 0; // Intervall从0开始
foreach(输入中的var行)
{
if(string.IsNullOrWhiteSpace(line))continue; //删除空行

x + = UsedInterval; //添加intervall
double y = double.Parse(line.Replace(,,。),CultureInfo.InvariantCulture);

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

我希望你能帮助我解决我的问题,并感谢你的关注。 / p>

当我做一个256位的ZeroPadding FFT时,上面的数据在一个周期内,间隔0,00001953125 这是我的输出 复杂输入与 Sin(coeff * i)以及零 - 纯真实输入的虚数部分,您将得到正确的输出,没有接近零的大峰值 - 这是由于当前实数部分的线性分量。

some days ago I posted my FFT-Test tool with some problems, so this here is my "remake" poste with almost the same problems, but different code.

Problem: I´m trying to make a tool to read some sensor values and do the fft of these. I implemented a normal DFT, a FFT(Bluestein) and a FFT(Cooley-Tukey). It all works fine and the inverse transformation shows the input like it should. The only thing is when i plot the Frequenz/Amplitude graph i have a huge peak at the begining (not only at 0 like the DC (0Hz) component should be!).

For Bluestein and DFT I use the complet input and for Cooley-Tukey I calculate with ZeroPadding or just cut the input to a 2^n length.

For example here is my input, just a simple sinuswave. Here the output with DFT and Bluetstein(looks the same) and FFT(ZeroPadding) and one DFT zoomed in at 0Hz. And here is the ouput of the same input with octave, what I expected to get.

Here my Code fore 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;
            }
        }

And here Bluestein.

            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;
        }

And this is the input per periode with an interval of 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

I have no idear what this peak could be. The only thing i can say is, when my interval is bigger the peak at the beginning is also bigger: Interval 0,00001953125 ; 0,0001953125 ; 0,001953125

Here is how i create the input of the fft. Just my shown input and the intervall....

                    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));
                }

I hope you can help me with my problems and thanks for your attention.

When I do a 256bit ZeroPadding FFT with the upper Data over one periode and the interval 0,00001953125 this is my ouput

解决方案

You should fill real part of complex input with Sin(coeff*i) and imaginary part with zero - pure real input, and you'll get right output without large peak near zero - it is due to linear component of current real part.

这篇关于FFT(Bluestein和Cooley-Tukey)...在开始时意外的高峰的文章就介绍到这了,希望我们推荐的答案对大家有所帮助,也希望大家多多支持IT屋!

查看全文
登录 关闭
扫码关注1秒登录
发送“验证码”获取 | 15天全站免登陆