Matlab FFT和自制FFT [英] Matlab FFT and home brewed FFT

问题描述

我正在尝试验证在VS Matlab上应该用于同一项目的FFT算法. 关键是，有了我自己的C FFT函数，我总是得到在Matlab中评估的双面FFT频谱的正确部分(第二个部分)，而不是预期"的第一个部分.

I'm trying to verify an FFT algorithm I should use for a project VS the same thing on Matlab. The point is that with my own C FFT function I always get the right (the second one) part of the double sided FFT spectrum evaluated in Matlab and not the first one as "expected".

例如，如果我的第三个bin的格式为a + i * b，则Matlab FFT的第三个bin的格式为a-i * b. A和B值相同，但我总是得到Matlab的复共轭. 我知道在振幅和功率方面没有问题(因为有绝对值)，但我想知道在相位方面我是否会总是读取错误的角度.

For instance if my third bin is in the form a+i*b the third bin of Matlab's FFT is a-i*b. A and b values are the same but i always get the complex conjugate of Matlab's. I know that in terms of amplitudes and power there's no trouble (cause abs value) but I wonder if in terms of phases I'm going to read always wrong angles.

我对Matlab不太了解(并且我在网上没有找到有用的信息)Matlab FFT是否可能首先以负频率返回FFT频谱然后再以正频率返回...或者我是否必须修复FFT算法...或者是否还可以，因为无论FFT的哪个部分，相位都是不变的，所以我们选择作为单边频谱(但是我对最后一个选择表示怀疑).

Im not so skilled in Matlab to know (and I have not found useful infos on the web) if Matlab FFT maybe returns the FFT spectre with negative frequencies first and then positive... or if I have to fix my FFT algorithm... or if it is all ok because phases are the unchanged regardless wich part of FFT we choose as single side spectrum (but i doubt about this last option).

示例:

如果S是具有N = 512个样本的样本阵列，则Matlab中的Y = fft(S)将FFT返回为(阵列前半部分的虚部的符号是随机的，只是为了显示复共轭第二部分的差异):

If S is the sample array with N=512 samples, Y = fft(S) in Matlab return the FFT as (the sign of the imaginary part in the first half of the array are random, just to show the complex conjugate difference for the second part):

1   A1    +  i*B1     (DC, B1 is always zero)
2   A2    +  i*B2
3   A3    -  i*B3
4   A4    +  i*B4
5   A5    +  i*B5
...
253 A253  -  i*B253
254 A254  +  i*B254
255 A255  +  i*B255
256 A256  +  i*B256
257 A257  -  i*B257   (Nyquyst, B257 is always zero)
258 A256  -  i*B256
259 A255  -  i*B255
260 A254  -  i*B254
261 A253  +  i*B253
...
509 A5    -  i*B5
510 A4    -  i*B4
511 A3    +  i*B3
512 A2    -  i*B2

我的FFT实现在Y数组中仅返回256个值(没关系)，

My FFT implementation returns only 256 values (and that's ok) in the the Y array as:

1   1    A1    +  i*B1     (A1 is the DC, B1 is Nyquist, both are pure Real numbers)
2   512  A2    -  i*B2
3   511  A3    +  i*B3
4   510  A4    -  i*B4
5   509  A5    +  i*B5
...
253 261  A253  +  i*B253
254 260  A254  -  i*B254
255 259  A255  -  i*B255
256 258  A256  -  i*B256

第一列是我的Y数组的正确索引，第二列只是Matlab FFT实现中相对行的引用.

Where the first column is the proper index of my Y array and the second is just the reference of the relative row in the Matlab FFT implementation.

如您所见，我的FFT实现(DC分开)返回FFT，就像 Matlab FFT的后半部分一样(相反顺序).

总结:即使我按照建议使用fftshift，看来我的实现总是返回在Matlab FFT中应该被认为是频谱的负部分. 错误在哪里?

To summarize: even if I use fftshift as suggested, it seems that my implementation always return what in the Matlab FFT should be considered the negative part of the spectrum. Where is the error???

这是我使用的代码:

注1:此处未声明FFT数组，并且在函数内部对其进行了更改. 最初包含N个样本(真实值)，最后 包含单边FFT频谱的N/2 +1个单元

Note 1: the FFT array is not declared here and it is changed inside the function. Initially it holds the N samples (real values) and at the end it contains the N/2 +1 bins of the single sided FFT spectrum.

注2:N/2 + 1个二进制位存储在N/2个元素中，仅是因为DC分量始终是实数(并且存储在FFT [0]中)，同时也包括Nyquyst(并且它存储在FFT中) [1])，除了所有其他偶数元素K都拥有一个实数，而烤箱元素K + 1则拥有虚部，这一例外.

Note 2: the N/2+1 bins are stored in N/2 elements only because the DC component is always real (and it is stored in FFT[0]) and also the Nyquyst (and it is stored in FFT[1]), this exception apart all the other even elements K holds a real number and the oven elements K+1 holds the imaginary part.

void Fft::FastFourierTransform( bool inverseFft ) {
   double twr, twi, twpr, twpi, twtemp, ttheta;
   int i, i1, i2, i3, i4, c1, c2;
   double h1r, h1i, h2r, h2i, wrs, wis;
   int nn, ii, jj, n, mmax, m, j, istep, isign;
   double wtemp, wr, wpr, wpi, wi;    
   double theta, tempr, tempi;

   // NS is the number of samples and it must be a power of two
   if( NS == 1 )
      return;

   if( !inverseFft ) {
      ttheta = 2.0 * PI / NS;
      c1 = 0.5;
      c2 = -0.5;
   }
   else {
      ttheta = 2.0 * PI / NS;
      c1 = 0.5;
      c2 = 0.5;
      ttheta = -ttheta;
      twpr = -2.0 * Pow( Sin( 0.5 * ttheta ), 2 );
      twpi = Sin(ttheta);
      twr = 1.0+twpr;
      twi = twpi;
      for( i = 2; i <= NS/4+1; i++ ) {
         i1 = i+i-2;
         i2 = i1+1;
         i3 = NS+1-i2;
         i4 = i3+1;
         wrs = twr;
         wis = twi;
         h1r = c1*(FFT[i1]+FFT[i3]);
         h1i = c1*(FFT[i2]-FFT[i4]);
         h2r = -c2*(FFT[i2]+FFT[i4]);
         h2i = c2*(FFT[i1]-FFT[i3]);
         FFT[i1] = h1r+wrs*h2r-wis*h2i;
         FFT[i2] = h1i+wrs*h2i+wis*h2r;
         FFT[i3] = h1r-wrs*h2r+wis*h2i;
         FFT[i4] = -h1i+wrs*h2i+wis*h2r;
         twtemp = twr;
         twr = twr*twpr-twi*twpi+twr;
         twi = twi*twpr+twtemp*twpi+twi;
      }
      h1r = FFT[0];
      FFT[0] = c1*(h1r+FFT[1]);
      FFT[1] = c1*(h1r-FFT[1]);
   }
   if( inverseFft )
      isign = -1;
   else
      isign = 1;
   n = NS;
   nn = NS/2;
   j = 1;
   for(ii = 1; ii <= nn; ii++) {
      i = 2*ii-1;
      if( j>i ) {
         tempr = FFT[j-1];
         tempi = FFT[j];
         FFT[j-1] = FFT[i-1];
         FFT[j] = FFT[i];
         FFT[i-1] = tempr;
         FFT[i] = tempi;
      }
      m = n/2;
      while( m>=2 && j>m ) {
         j = j-m;
         m = m/2;
      }
      j = j+m;
   }
   mmax = 2;
   while(n>mmax) {
      istep = 2*mmax;
      theta = 2.0 * PI /(isign*mmax);
      wpr = -2.0 * Pow( Sin( 0.5 * theta ), 2 );
      wpi = Sin(theta);
      wr = 1.0;
      wi = 0.0;
      for(ii = 1; ii <= mmax/2; ii++) {
         m = 2*ii-1;
         for(jj = 0; jj <= (n-m)/istep; jj++) {
            i = m+jj*istep;
            j = i+mmax;
            tempr = wr*FFT[j-1]-wi*FFT[j];
            tempi = wr*FFT[j]+wi*FFT[j-1];
            FFT[j-1] = FFT[i-1]-tempr;
            FFT[j] = FFT[i]-tempi;
            FFT[i-1] = FFT[i-1]+tempr;
            FFT[i] = FFT[i]+tempi;
         }
         wtemp = wr;
         wr = wr*wpr-wi*wpi+wr;
         wi = wi*wpr+wtemp*wpi+wi;
      }
      mmax = istep;
   }
   if( inverseFft )
      for(i = 1; i <= 2*nn; i++)
         FFT[i-1] = FFT[i-1]/nn;
   if( !inverseFft ) {
      twpr = -2.0 * Pow( Sin( 0.5 * ttheta ), 2 );
      twpi = Sin(ttheta);
      twr = 1.0+twpr;
      twi = twpi;
      for(i = 2; i <= NS/4+1; i++) {
         i1 = i+i-2;
         i2 = i1+1;
         i3 = NS+1-i2;
         i4 = i3+1;
         wrs = twr;
         wis = twi;
         h1r = c1*(FFT[i1]+FFT[i3]);
         h1i = c1*(FFT[i2]-FFT[i4]);
         h2r = -c2*(FFT[i2]+FFT[i4]);
         h2i = c2*(FFT[i1]-FFT[i3]);
         FFT[i1] = h1r+wrs*h2r-wis*h2i;
         FFT[i2] = h1i+wrs*h2i+wis*h2r;
         FFT[i3] = h1r-wrs*h2r+wis*h2i;
         FFT[i4] = -h1i+wrs*h2i+wis*h2r;
         twtemp = twr;
         twr = twr*twpr-twi*twpi+twr;
         twi = twi*twpr+twtemp*twpi+twi;
      }
      h1r = FFT[0];
      FFT[0] = h1r+FFT[1];  // DC
      FFT[1] = h1r-FFT[1];  // FS/2 (NYQUIST)
   }            
   return;
}

推荐答案

在matlab中，尝试使用fftshift(fft(...)).在调用FFT后，Matlab不会自动移动频谱，这就是为什么他们实现了fftshift()函数的原因.

In matlab try using fftshift(fft(...)). Matlab doesn't automatically shift the spectrum after the FFT is called which is why they implemented the fftshift() function.

这篇关于Matlab FFT和自制FFT的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！