Matlab FFT与国产FFT

Question

我试图验证一个FFT算法，我应该使用一个项目和相同的东西在Matlab上。问题是，用我自己的C FFT函数，我总是得到在Matlab中评估的双面FFT谱的正确部分(第二部分)，而不是像“预期”那样得到第一个部分。

例如，如果我的第三个宾子是a+i*b，那么Matlab的FFT的第三个值是a-i*b。a和b值是相同的，但是我总是得到Matlab的复共轭。我知道在振幅和功率方面没有问题(引起abs值)，但是我不知道从相位上看，我是否总是读取错误的角度。

我不太擅长在Matlab知道(我还没有在网上找到有用的信息)，如果Matlab FFT可能返回的负频率的FFT频谱，然后正…或者如果我要修正FFT算法..。或者，如果一切正常，因为相位是不变的，无论FFT的哪个部分，我们选择作为单边频谱(但我怀疑这最后的选择)。

示例：

如果S是带有N=512样本的样本阵列，则在Matlab中Y= fft(S)返回FFT作为(阵列前半部分虚部的符号是随机的，只是为了表示第二部分的复共轭差)：

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个值(这很好)，如下所示：

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实现(DC分离)返回FFT，就像Matlab的FFT的后半部分(以相反的顺序)。

总之:即使我按建议使用fftshift，我的实现似乎总是返回Matlab中应该被认为是频谱中的负面部分。错误在哪里？

，这是我使用的代码：

注1:这里没有声明FFT数组，而是在函数中对其进行了更改。它最初保存N个样本(实值)，最后包含单边FFT谱的N/2 +1桶。

注2: N/2+1桶被存储在N/2元素中，这是因为DC分量总是真实的(并且存储在FFT中)和Nyquyst (并且存储在FFT1中)，除了所有其他偶数元素K保持一个实数，而K+1保存虚部。

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

Answer 1

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

Answer 2

它只是一个matlab格式的东西。基本上，matlab按如下顺序排列傅里叶变换。

DC, (DC-1), .... (Nyquist-1), -Nyquist, -Nyquist+1, ..., DC-1

假设你有一个8点序列:1 2 3 1 4 5 1 3

在您的信号处理类中，您的教授可能会根据笛卡尔系统(x轴为负->正)绘制傅里叶谱；因此，您的DC应该位于x轴上的0(在您的快速傅立叶变换序列中的第四个位置，假设位置指数为0)。

在matlab中，DC是fft序列中的第一个元素，因此您需要用fftshit()来交换fft序列的前半部分和后半部分，以便DC位于第4位置(位置为0基索引)。

我在这里附加了一个图表，这样你就可以看到：

​

​

其中a是原始的8点序列，FT(a)是a的傅里叶变换.

matlab代码在这里：

a = [1 2 3 1 4 5 1 3];
A = fft(a);

N = length(a);
x = -N/2:N/2-1;

figure
subplot(3,1,1), stem(x, a,'o'); title('a'); xlabel('time')
subplot(3,1,2), stem(x, fftshift(abs(A),2),'o'); title('FT(a) in signal processing'); xlabel('frequency')
subplot(3,1,3), stem(x, abs(A),'o'); title('FT(a) in matlab'); xlabel('frequency')