2016-09-21 123 views
-4

我目前正处于实习阶段,我的老板希望我能够在本周末用java编写准备好FFT。了解FFT输出

现在我已经得到了来自普林斯顿大学的FFT代码:http://introcs.cs.princeton.edu/java/97data/FFT.java

我实现和扩展这个代码在我的项目与我现在能先读一个信号的二进制输入,然后在这个采样值上处理FFT,然后提供幅度。

现在我来找我的问题。 我输入了我生成的以下值,并使用final double d = Math.sin(i);并循环了8次(这仅用于测试目的,下周我将不得不输入实际值)。

0.0 
0.8414709848078965 
0.9092974268256817 
0.1411200080598672 
-0.7568024953079282 
-0.9589242746631385 
-0.27941549819892586 
0.6569865987187891 

因此这些值来自纯正弦波(我不知道正确的英文单词,但与纯正弦波我的意思是可以精确到一个频率,例如50赫兹正弦)。

的输出现在是

0.553732750242242 
2.3946469565385193 - 2.0970118573701813i 
-1.386684423934684 + 0.9155598966338983i 
-0.8810419659226628 + 0.28041399267903344i 
-0.8075738836045867 
-0.8810419659226628 - 0.28041399267903366i 
-1.386684423934684 - 0.9155598966338983i 
2.394646956538519 + 2.0970118573701817i 

和输出

0.553732750242242 
3.183047718211326 
1.6616689248786416 
0.9245901540720989 
0.8075738836045867 
0.924590154072099 
1.6616689248786416 
3.183047718211326 

现在的大小实际上我预期的输出值在每个频率样本为0点,直到我到达频率的纯正弦由输出应该> 0(例如50Hz)控制。至少这是我的老板在给我这个任务时所期望的。

总结: 所以这就是我正在努力。我读过另一个线索,询问有关类似的问题,但对我来说还有一些未解决的问题。我该如何处理给定的输出数据?我如何找到最频繁的频率?

我真的可能需要一些帮助或解释我的思维错误。

感谢收听...

回答

1

计算应用简单的窗函数后512点的傅立叶变换:

w(i)= ((float)i-(float)(n-1f)/2f) 

它在I = 25(结果阵列上最大大小)给出峰值。还加入

输入用的详细信息,例如正弦波发生器(50赫兹)和采样率(1kHz的或每个样品0.001秒)的频率和添加2PI常数:

初始化现在看起来像这样(作为罪(2xPIxFxi)表示法):

for (int i = 0; i < n; i++) 
{ 
     v1[i] = ((float)i-(float)(n-1f)/2f)* 
      (float)Math.Sin(0.001f* 50f*2f*Math.PI*(float)i); 
          ^ ^^    ^
           |  | |     | 
           |  F |     | 
          sampling  2xPi constant sample bin 
          rate(simulation) 
} 

,并导致峰的样子:

v2[22] 2145,21852033773 
v2[23] 3283,36245333956 
v2[24] 6368,06249969329 
v2[25] 28160,6579468591 <-- peak 
v2[26] 23231,0481898687 
v2[27] 1503,8455705291 
v2[28] 1708,68502071037 

所以我们必须

25 

现在这些步骤在频率空间中,输入频率为1kHz,因此您可以感知最大500 Hz的谐波信号(采样率的一半)。

25*500 = 12500 

也导致范围为0到N/2与另一半镜像,以及将可察觉的频率范围(500)的结果的范围(256对于N = 512)给出的

48.83 Hz 

大部分误差一定是在开始时使用的窗口函数,但v2 [26]的值比v2更高[24],因此实际采样值更接近v2 [26],这些点的平滑图应显示50 Hz

忽略结果数组的第一个元素,因为它关于恒定信号电平或无限远波长或零频率。

下面是DFT计算的代码只是如果FFT将返回正确的结果是肯定的:

//a---->b Fourier Transformation brute-force 
__kernel void dft(__global float *aRe, 
        __global float *aIm, 
        __global float *bRe, 
        __global float *bIm) 
{ 
    int id=get_global_id(0); // thread id 
    int s=get_global_size(0); // total threads = 512 
    double cRe=0.0f; 
    double cIm=0.0f; 
    double fid=(double)id; 
    double fmpi2n=(-2.0*M_PI)*fid/(double)s; 
    for(int i=0;i<s;i++) 
    { 
      double fi=(float)i; 
      double re=cos(fmpi2n*fi); 
      double im=sin(fmpi2n*fi); 

      cRe+=aRe[i]*re-aIm[i]*im; 
      cIm+=aRe[i]*im+aIm[i]*re; 
    } 

    bRe[id]=cRe; 
    bIm[id]=cIm; 
} 

,并予以肯定,对逆变换的测试结果,以检查是否原始输入信号再次实现:

// a--->b inverse Fourier Transformation brute force 
__kernel void idft(__global float *aRe, 
        __global float *aIm, 
        __global float *bRe, 
        __global float *bIm) 
{ 
    int id=get_global_id(0); // thread id 
    int s=get_global_size(0); // total threads = 512 
    double cRe=0.0f; 
    double cIm=0.0f; 
    for(int i=0;i<s;i++) 
    { 
     double re=cos(2.0*M_PI*((double)id)*((double)i)/(double)s); 
     double im=sin(2.0*M_PI*((double)id)*((double)i)/(double)s); 

     cRe+=aRe[i]*re-aIm[i]*im; 
     cIm+=aRe[i]*im+aIm[i]*re; 
    } 
    cRe/=(double)s; 
    cIm/=(double)s; 
    bRe[id]=cRe; 
    bIm[id]=cIm; 
} 

我知道在快速机器上运行速度慢的代码是不好的,但是这看起来好简单多了,可以扩展到许多内核(对于包含阵列副本的320core-gpu是2.4ms)。

+0

你好侯赛因,非常感谢,这非常有帮助! – Winterwurst

+0

@Winterwurst如果峰值不在您预见的位置,您可以尝试插值以精确调整两个整数之间的峰值点。 –