fft算法.txt

我收集的fft算法

    // 中间变量
    int        bfsize,p  ;
    
    // 角度
    float    angle;
    
    complex<float> *W,*X1,*X2,*X;
    
    // 计算付立叶变换点数  
    count = 1 << r;
    
    // 分配运算所需存储器
    W  = new complex<float>[count / 2];
    X1 = new complex<float>[count];
    X2 = new complex<float>[count];
    
    // 计算加权系数
    for(i = 0; i < count / 2; i++)
    {
        angle = -i * PI * 2 / count;
        W[i] = complex<float> (cos(angle), sin(angle));
    }
    
    // 将时域点写入X1
    memcpy(X1, TD, sizeof(complex<float>) * count);
    
    // 采用蝶形算法进行快速付立叶变换
    for(k = 0; k < r; k++)
    {
        for(j = 0; j < 1 << k; j++)
        {
            bfsize = 1 << (r-k);
            for(i = 0; i < bfsize / 2; i++)
            {
                p = j * bfsize;
                X2[i + p] = X1[i + p] + X1[i + p + bfsize / 2];
                X2[i + p + bfsize / 2] = (X1[i + p] - X1[i + p + bfsize / 2]) * W[i * (1<<k)];
            }
        }
        X  = X1;
        X1 = X2;
        X2 = X;
    }
    
    // 重新排序
    for(j = 0; j < count; j++)
    {
        p = 0;
        for(i = 0; i < r; i++)
        {
            if (j&(1<<i))
            {
                p+=1<<(r-i-1);
            }
        }
        FD[j]=X1[p];
    }
    
    // 释放内存
    delete W;
    delete X1;
    delete X2;
} 



cos(x)的定点实现
用泰勒公式的话，要求x^2,x^4,x^6,......,这样运算量很大（要求高精度的话）.有没有别的办法？
 

yangrui
我是一个兵  第4楼 回复于2004-9-1 17:47:44  
--------------------------------------------------------------------------------

这是我用定点做的，希望大家提意见！！！！

//Fix Point cosine(x) function  in/out Q15 
Word16 Cosine(Word16 x)   
{
    /* table of cos(x) in Q15 */
    Word16 table[65] = {
          32767,  32729,  32610,  32413,  32138,  31786,  31357,  30853,
          30274,  29622,  28899,  28106,  27246,  26320,  25330,  24279,
          23170,  22006,  20788,  19520,  18205,  16846,  15447,  14010,
          12540,  11039,   9512,   7962,   6393,   4808,   3212,   1608,
              0,  -1608,  -3212,  -4808,  -6393,  -7962,  -9512, -11039,
         -12540, -14010, -15447, -16846, -18205, -19520, -20788, -22006,
         -23170, -24279, -25330, -26320, -27246, -28106, -28899, -29622,
         -30274, -30853, -31357, -31786, -32138, -32413, -32610, -32729,
         -32768L };

    /* slop in Q15 */
    Word16 slope_cos[64] = {
        -774,   -2424,  -4013,  -5602,  -7170,  -8739,  -10267, -11795, 
        -13282, -14728, -16154, -17519, -18864, -20168, -21410, -22592, 
        -23712, -24812, -25831, -26788, -27685, -28500, -29274, -29946, 
        -30578, -31107, -31576, -31963, -32289, -32513, -32676, -32757, 
        -32757, -32676, -32513, -32289, -31963, -31576, -31107, -30578, 
        -29946, -29274, -28500, -27685, -26788, -25831, -24812, -23712, 
        -22592, -21410, -20168, -18864, -17519, -16154, -14728, -13282, 
        -11795, -10267, -8739,  -7170,   -5602, -4013,  -2424,  -774
    }; 

    int inp;

    Word16 step;
    Word32 L_tmp;
    Word16 val;

    inp = 0;
    step = 0;
    while((x > step) && (inp < 64))
    {
        step += 1608;    /* 1608 is PI/64 in Q15 */
        inp++;
    }

    if(inp >= 64)
        inp = 63;

    L_tmp  = L_mult(sub(x, step) , slope_cos[inp] );   /* L_tmp Q31 */
         
    val = add(table[inp], extract_h(L_tmp));           /* return Q15 */

    return val;
}
 

要求高精度，运算量肯定是很大的，
泰勒公式肯定可以做到，就看你怎么去实现了，
比如，x^6 = x^4 * x^2之类的小技巧等 




正弦波的傅立叶变换应该是脉冲是没错的。
你如果对一个正铉波在一个周期内取样1024个点，假设正铉波的频率是1Hz，那么取样频率就是1024Hz，fft后，得到1024个复数点。计算出这些复数点的幅度值，那么这1024个点就对应相应的频率。
因为fft是对称的，所以结果从0-511就对应了1-512Hz,频谱分辨率是1Hz,在1Hz的地方，肯定有个脉冲。
可以在Matlab运行以下程序：
t=linspace(1,1024,1024);
y=sin(t.*(2*pi/1024));
psd=abs(fft(y));
plot(linspace(0,19,20),psd(1:20)); 
xlabel('Frequency Hz'); 



我这儿也有一个c版的源程序
[代码]FFT 快速傅立叶变换 (c语言)

算法参见徐世良《计算机常用算法》第二版 第294页

void KFFT(
BNU::vector<std::complex<double> > &x,
BNU::vector<std::complex<double> > &y,
int inverse)
{
  int n = x.size();
  int k = log((double)n)/log(2.);
  int t ,m ,s , i, j, NV, r;
  double angle;
  std::complex<double> V,PODD;
  for (t = 0 ; t <= n-1 ; t++) {
    m = t;
    s = 0;
    for (i = 0 ; i <=k-1 ; i++) {
      j = int(m/2);
      s = 2*s + m - 2*j;
      m = j;
    }                                         
  if (&x!=&y) {
      y[t] = x[s];
    }else {
       if (s>t) {
           V=x[t];
           y[t]=x[s];
           y[s]=V;
       }
    }
  }
  for (t = 0 ; t <= n-2 ; t+=2) {
    V = y[t];
    y[t] = V + y[t+1];
    y[t+1] = V - y[t+1];
  }
  m = n/2;
  NV = 2;
  for (r = k-2 ; r >=0 ; r-- ) {
    m = m/2;
    NV = 2*NV;
    for (t = 0 ; t <= (m-1)*NV ; t=t+NV ){
      for (j = 0 ; j <= (NV/2)-1 ; j++) {             
  angle = -(2.0*M_PI*inverse*j)/NV;
        PODD = complex<double>(cos(angle),sin(angle))*y[t+j+(NV/2)];
        y[t+j+(NV/2)] = y[t+j] - PODD;
        y[t+j] = y[t+j] + PODD;
      }
    }
  };
for (i = 0 ; i < n; i++) {
   if (inverse==-1)
        y = y/double(n);
}
} 




____________________________________________________________________



Dim num1, num2, num3
Dim a(), a1(), e() As Single
Dim alf, xgm(), p(), k1(), k2, k(), p1(), p2(), c(), c1, k3(), k4(), cc() As Single
Private Sub Command1_Click()
CommonDialog1.ShowOpen
Open CommonDialog1.FileName For Input As #1
num3 = Val(Text1.Text) + 2
num2 = Val(Text2.Text)
ReDim a(num3, num2)
For j = 1 To num3
For i = 1 To num2
Input #1, a(j, i)
Next i
Next j
Close #1

End Sub

Private Sub Command2_Click()
CommonDialog1.ShowSave
Open CommonDialog1.FileName For Append As #2

For i = 1 To num2
Print #2, cc(i, 1), cc(i, 2), cc(i, 3)
Next i

Close #2
Shell "c:\windows\notepad " & CommonDialog1.FileName
End Sub

Private Sub Command3_Click()

num1 = Val(Text1.Text) + 1
num2 = Val(Text2.Text)
num3 = Val(Text1.Text) + 2
ReDim a1(num2), e(num1, num2) As Single
ReDim xgm(num1), p(num1, num1), k1(num1), k(num1), p1(num1, num1), p2(num1, num1), c(num1), k3(num1), k4(num1), cc(num2, num1) As Single





c(1) = c(2) = c(3) = 0
r = 0.000001
alf = 10
For j = 1 To num1
For i = 1 To num2
e(j, i) = a(j, i)
a1(i) = a(num3, i)
Next i
Next j
For x = 1 To num1
xgm(x) = alf * Sqr(r / e(x, 1))
p(x, x) = xgm(x) ^ 2
Next x

For q = 1 To num2
For i = 1 To num1
k3(i) = k4(i) = 0
Next i
k2 = 0
For i = 1 To num1
For j = 1 To num1
p1(i, j) = p2(i, j) = 0
Next j
Next i
c1 = 0

For m = 1 To num1
k3(m) = 0
For n = 1 To num1
k3(m) = k3(m) + e(n, q) * p(n, m)
Next n
Next m

For m = 1 To num1
k4(m) = 0
k4(m) = k4(m) + 1000 * k3(m) * e(m, q)
Next m

For i = 1 To num1
k2 = k2 + k4(i)
Next i

For v = 1 To num1
k(v) = 0
k(v) = k3(v) * (k2 / 1000 + r) ^ -1
Next v

For i = 1 To num1
For j = 1 To 3
p1(i, j) = k(i) * e(j, q)
Next j
Next i

For l = 1 To num1
For i = 1 To num1
p2(l, i) = 0
For j = 1 To num1

p2(l, i) = 1000 * p1(l, j) * p(j, i) + p2(l, i)
Next j
Next i
Next l

For i = 1 To num1
For j = 1 To num1
p(i, j) = p(i, j) - p2(i, j) / 1000
Next j
Next i

For i = 1 To num1
c1 = c1 + e(i, q) * c(i)
Next i
For i = 1 To num1
c(i) = c(i) + k(i) * (a1(q) - c1)

Next i


For i = 1 To num1

cc(q, i) = c(i)
Next i
Next q

End Sub

Private Sub Command4_Click()
End
End Sub

上面是用vb编写的卡尔曼滤波算法程序（转贴），应该不难看懂 






由于需要，急需积的算法，谁能帮帮我？  

happyyangxu
纯净水  第1楼 回复于2004-6-11 21:18:05    
--------------------------------------------------------------------------------
 本回复被接受作为正确答案


我有个给你吧！
float f1(float x)
{
　float f;
　f = 1 + x*x;
　return f;
}

float f2(float x)
{
　float f;
　f = 1 + x + x*x + x*x*x;
　return f;
}

float f3(float x)
{
　float f;
　f = x / (1 + x*x);
　return f;
}



float integral(float (*fun)(float), float a, float b)  
{
　float s, h, y;
　int n, i;
　
　s = ( (*fun)(a) + (*fun)(b) ) /2.0;
　n = 100;
　h = (b-a)/n;
　for(i=1; i<n; i++)
　　　s = s + (*fun)(a+i*h);
　y = s * h;
　return y;
} 
void main()
{
　float y1, y2, y3;
　y1 = integral(f1, 0.0, 1.0);
　y2 = integral(f2, 0.0, 2.0);
　y3 = integral(f3, 0.0, 3.5);
　printf("y1=%6.2fny2=%6.2fny3=%6.2fn", y1,y2,y3); 
}