表教高手解答一个关于Cfft问题

只看该作者 · 2009-12-28 11:52

由于这儿传图片太不方便了，我把它作成了WORD文档，希望高手看一下。谢谢。

只看该作者 · 2009-12-29 08:42

自己顶一下，难道没人指教一下吗？

只看该作者 · 2010-1-3 01:56

FFT 有两种实现方法：
1. 时域抽取要先bit reverse，然后butterfly运算。
2. 频域抽取时先做bufferfly运算，然后再bit reverse。

你的FFT应该是频域抽取。

建议你看一下discrete-time signal processing second edition, Oppenheim

只看该作者 · 2010-1-4 17:34

3# aresc 谢谢你的回复。我用的函数都是DSPLIB库中的函数。所以不知道它是频域抽取还是时域抽取。但是我还有一个地方不明白：在TI的文档里说的是CFFT函数的输入数据应该是复数啊。也就是需要两个存储空间才能表示一个数。但我用AD采样采出来的数应该都是实数啊。为什么直接赋给CFFT函数也能计算出来结果呢？

只看该作者 · 2010-1-4 22:32

实数只是复数的子集而已，虚部为0的复数不就是实数。

你应该好好看一下CFFT的输入是什么样的格式，c = a + bj，实部和虚部是怎样传进去的。

你看不到CFFT函数的具体实现代码吗？没有调用说明吗？

只看该作者 · 2010-1-5 10:20

5# aresc
再次感谢您的回复！TI文档是这样说的：
Function void cfft (DATA *x, ushort nx, type);
(defined in cfft.asm)
Arguments
x [2*nx] Pointer to input vector containing nx complex elements (2*nx
real elements) in normal order. On output, vector contains
the nx complex elements of the FFT(x) in bit-reversed order.
Complex numbers are stored in interleaved Re-Im format.
nx Number of complex elements in vector x. Must be between
16 and 1024.
type FFT type selector. Types supported:
If type = SCALE, scaled version selected
If type = NOSCALE, non-scaled version selected
我用的是5502，定点的DSP，数组X如果是复数的话，应该这样的一个存储结构：
Xr(N)
Xi(N)
Xr(N-1)
Xi(N-1)
...
一个实部，一个虚部的结构。我迷惑的是：我用AD采样采的都是实数啊，直接就赋值给了X[]，虚部也没有置0啊，为什么得出来的结果也是对的？是不是这句话有问题：
x [2*nx] Pointer to input vector containing nx complex elements (2*nx
real elements) in normal order.我理解的是：X数组既可以传入NX个复数，也可以传入2*NX个实数，它们都可以进行FFT的变换，不知道我这样理解对不对？

x

只看该作者 · 2010-1-5 10:23

5# aresc
这是CCS给的例子：
/* Test for cfft: rand */
#define NX 64
#define FNAME "t4"
#define MAXERROR 16
#pragma DATA_SECTION (x,".input")
DATA x[2*NX] ={
13971, 9823, /* 0 */
8918, 11219, /* 1 */
-1919, -1066, /* 2 */
3344, 2704, /* 3 */
14508, -18688, /* 4 */
-21739, -2429, /* 5 */
5072, 10645, /* 6 */
-22178, 19372, /* 7 */
-4174, -11190, /* 8 */
-9617, -12094, /* 9 */
18461, 18029, /* 10 */
-23917, -13199, /* 11 */
13213, 15034, /* 12 */
23219, 20139, /* 13 */
24167, -13222, /* 14 */
14244, -12856, /* 15 */
-3025, -22204, /* 16 */
-84, -20792, /* 17 */
-14106, 6944, /* 18 */
7076, -15244, /* 19 */
-8875, 16957, /* 20 */
22689, -16082, /* 21 */
11176, -16235, /* 22 */
-4342, 24375, /* 23 */
12060, -2970, /* 24 */
-11444, -7888, /* 25 */
-2963, -9162, /* 26 */
21371, -6654, /* 27 */
9040, -5265, /* 28 */
-14175, 4513, /* 29 */
16729, -18752, /* 30 */
6350, -22777, /* 31 */
-18060, -2042, /* 32 */
-14443, 18239, /* 33 */
5286, 21413, /* 34 */
6405, -11616, /* 35 */
-6388, -16752, /* 36 */
3705, 18386, /* 37 */
-2396, -12927, /* 38 */
-22493, 7191, /* 39 */
-23317, 23023, /* 40 */
-9238, 8133, /* 41 */
-24023, 18264, /* 42 */
-5722, -24168, /* 43 */
9030, -17901, /* 44 */
-20079, 15718, /* 45 */
-22915, -3444, /* 46 */
5542, 19248, /* 47 */
5352, 11584, /* 48 */
-23880, 9237, /* 49 */
-23851, -7589, /* 50 */
-15284, -16470, /* 51 */
4286, -16983, /* 52 */
-21818, -15233, /* 53 */
-6531, -3825, /* 54 */
6482, 17554, /* 55 */
10732, -481, /* 56 */
9501, 15579, /* 57 */
-20511, -1935, /* 58 */
-2251, -2104, /* 59 */
-2869, -2432, /* 60 */
-7237, -4329, /* 61 */
-17082, 19804, /* 62 */
8661, -24382, /* 63 */
};
#pragma DATA_SECTION (rtest,".test")
DATA rtest[2*NX] ={
-2068, -785,  /* 0 */
-987, -3092,  /* 1 */
2094, 1909,  /* 2 */
-524, 1146,  /* 3 */
521, -491,  /* 4 */
-710, -905,  /* 5 */
1666, -2966,  /* 6 */
3081, -1402,  /* 7 */
-711, 3381,  /* 8 */
-1190, 21,  /* 9 */
455, 9,  /* 10 */
-177, 1146,  /* 11 */
2740, -657,  /* 12 */
-600, -30,  /* 13 */
2346, 212,  /* 14 */
1479, -4191,  /* 15 */
2847, 337,  /* 16 */
-218, 1229,  /* 17 */
-1663, -111,  /* 18 */
2344, 507,  /* 19 */
1728, 2919,  /* 20 */
1526, -74,  /* 21 */
142, -1693,  /* 22 */
1800, 282,  /* 23 */
-2999, -1415,  /* 24 */
310, -121,  /* 25 */
3505, 2798,  /* 26 */
1129, -233,  /* 27 */
-544, 815,  /* 28 */
-747, 3438,  /* 29 */
-212, 1371,  /* 30 */
-226, 1165,  /* 31 */
1132, -264,  /* 32 */
1722, -95,  /* 33 */
262, -1499,  /* 34 */
-3127, -1644,  /* 35 */
-1743, 925,  /* 36 */
229, 387,  /* 37 */
-2583, 2491,  /* 38 */
-1533, 875,  /* 39 */
2099, 475,  /* 40 */
-3592, -1913,  /* 41 */
-1162, 1295,  /* 42 */
625, 1462,  /* 43 */
-1877, -518,  /* 44 */
3158, 2511,  /* 45 */
-1412, -3200,  /* 46 */
965, -1061,  /* 47 */
-320, -1820,  /* 48 */
3218, 676,  /* 49 */
-1641, 942,  /* 50 */
348, -1325,  /* 51 */
-1880, -2501,  /* 52 */
-1410, 1972,  /* 53 */
-1786, 680,  /* 54 */
5987, -331,  /* 55 */
-791, 780,  /* 56 */
-1156, 1208,  /* 57 */
-745, 3239,  /* 58 */
2787, 368,  /* 59 */
1418, -1898,  /* 60 */
0, -161,  /* 61 */
-877, -884,  /* 62 */
1489, 4100,  /* 63 */
};
DATA error;
主程序的调用：
// compute
cfft(x,NX, SCALE);
// cfft(x,NX, NOSCALE);
cbrev(x,x,NX);

只看该作者 · 2010-1-5 11:05

TI的文档说的很清楚， void cfft (DATA *x, ushort nx, type)中x 是2 * nx个数据，其中一个复数占用相邻的2个数据(interleaved Re-Im format)，实部在前，虚部在后。而且也说了FFT的结果是bit reverse order的，所以是按频域抽取的。

由此可见你之前的那个结果不对，主要是进去的复数数据不对。你改一下，如果虚部全为0的话，应该和你的rfft结果一致。

只看该作者 · 2010-1-5 18:30

8# aresc
谢谢你一直的回复！按你说的，我把虚部置为0，重新进行CFFT的运算，结果确实和RFFT运算结果一样。但在调试过程中，有时会出现FFT变换后左右不对称的情况，是不是我调试的原因呢？下面是我的截图。

只看该作者 · 2010-1-5 20:28

你那两个图里都显示的是DDataBuffer里的数据，也就是已经经过bit reverse之后正常顺序的频域系数，re0,im0,re1,im1,re2,im2,...re1023,im1023，你的横轴是2048个点，两个图的纵轴分别是什么，是怎么计算那个纵坐标的值的？

只看该作者 · 2010-1-6 09:53

10# aresc
谢谢你的回复！在这个图中，我没有对DDataBuffer中的数据进行求模运算再显示，而是直接在CCS显示出来。图中纵轴的值就是DDataBuffer中数据的值，也可以理解为幅度。我通过把DDataBuffer中的数据进行取模之后再进行显示，两个波形就完全对称了。这是我的程序：
         for(i=0;i<1024;i++)
         {
         DDataBuffer[m]=DataBuffer;//实部为采样数据，虚部为0
         m+=2;
         }
            cfft_SCALE(DDataBuffer,1024);
      cbrev(DDataBuffer,DDataBuffer,1024);
      m=0;
      for(i=0;i<1024;i++)
      {

      real=DDataBuffer[m];
      m++;
      image=DDataBuffer[m];
      m++;
      DataBuffer=sqrt((long)real*(long)real+(long)image*(long)image);
      }

DDataBuffer中显示不对称是不是因为没有取模的原因？我知道要计算信号的幅度时是要取模的。可是我不明白它为什么不对称？下面的附件是我的截图。

只看该作者 · 2010-1-6 11:28

不取模也就是sqrt(real^2 + imag^2)你怎么显示2048个实部和虚部间隔的数据？只是简单的把实部、虚部显示在2048个时间点上，充其量你只能直观的从图形上看出实部、和虚部的分布情况，其他的没啥意义。

你程序中sqrt(real^2 + imag^2)的结果是不是总是取正的? 例如sqrt(4) = 2，但-2实际也是sqrt(4)的根，这样正好差一个符号，反对称的。

只看该作者 · 2010-1-6 11:40

12# aresc
谢谢您的指教！让我对DSPLIB库中的CFFT，RFFT都有了很深的了解。再次感谢你！如果你方便的话加我好友吧。我QQ：543433704，自己在研究DSP，没有老师的指导。很是吃力。专业课也学的不太好。呵呵。。。

只看该作者 · 2010-1-6 11:53

13# adofu2008
对了，我还有一个DSPLIB库中的IIR的问题不明白，也想请教你。
Cascaded IIR Direct Form II Using 4 Coefficients per Biquad
Function ushort oflag = iircas4 (DATA *x, DATA *h, DATA *r, DATA *dbuffer, ushort nbiq,
ushort nx)
(defined in iir4cas4.asm)
Arguments
x [nx] Pointer to input data vector of size nx
h[4*nbiq] Pointer to filter coefficient vector with the following
format:
h = a11 a21 b21 b11 ....a1i a2i b2i b1i
where i is the biquad index (a21 is the a2 coefficient of
biquad 1). Pole (recursive) coefficients = a. Zero
(non-recursive) coefficients = b
r[nx] Pointer to output data vector of size nx. r can be equal
dbuffer[2*nbiq] Pointer to address of delay line d.
Each biquad has 2 delay line elements separated by
nbiq locations in the following format:
d1(n-1), d2(n-1),..di(n-1) d1(n-2), d2(n-2)...di(n-2)
where i is the biquad index (d2(n-1) is the (n-1)th
delay element for biquad 2).
In the case of multiple-buffering schemes, this array
should be initialized to 0 for the first block only. Between
consecutive blocks, the delay buffer preserves
the previous r output elements needed.
Memory alignment: this is a circular buffer and must
start in a k-bit boundary(that is, the k LSBs of the
starting address must be zeros) where k = log2
(2*nbiq).
nbiq Number of biquads
nx Number of elements of input and output vectors
oflag Overflow flag.
If oflag = 1, a 32-bit overflow has occurred
If oflag = 0, a 32-bit overflow has not occurred
它说的每二阶四个系数是什么意思？是不是只有a1,a2,b1,b2?那我如何通过MATLAB得到这样的系数呢？我用matlab总是得到的是a0,a1,a2,b0,b1,b2?我对MATLAB不熟悉，只会利用FDATOOL设计滤波器。查了好多的资料也没有查出来。问合众达公司都一个月了也没有回信。

than x.

只看该作者 · 2010-1-6 13:24

biquad怎样变成4个系数：

         b0 + b1*z^-1 + b2*z^-2
H(z) = ------------------------                               (Eq 1)
         a0 + a1*z^-1 + a2*z^-2

         (b0/a0) + (b1/a0)*z^-1 + (b2/a0)*z^-2
H(z) = ---------------------------------------                (Eq 2)
            1 + (a1/a0)*z^-1 + (a2/a0)*z^-2

                           1 + (b1/b0)*z^-1 + (b2/b0)*z^-2
H(z) = (b0/a0) * ---------------------------------             (Eq 3)
                           1 + (a1/a0)*z^-1 + (a2/a0)*z^-2

推荐你看一本书 digital signal processing using matlab. 第6章有关于滤波器结构的。