FFT修改
FFT-dit的算法讨论
级的概念
M = l o g 2 N M=log_2N M=log2N,所以N点DFT可分成M级。
下图为8点FFT时间抽取算法信号流程图,方框分别为m=0,m=1,m=2;
码位倒置
交换后输出序列x(k)依照正序排列,但输入序列x(n)的次序不再是原来的自然次序,这是由于将x(n)按奇,偶分开产生的。
雷德算法: 对于自然顺序(二进制)我们是在低位加 1 得到下一位数,对于倒位序我们是在高位加 1 向低位进位。比如已知一个倒位序数是J求其下一个倒位序数,N位总数 ,把J与N/2比较若J<N/2则J的最高位为 0 ,把最高位置 1 ,就得到了J的下一个倒位序数;若J>=N/2则说明J的最高位为1 ,把最高位置0 ,比较次高位,若次高位为0 ,则把次高位置1,得到J的下一个倒位序,若次高位为1 , 则把次高位置0,以此类推…
参考:https://blog.csdn.net/corcplusplusorjava/article/details/17119567
Rader参考程序:
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
int main(void)
{
int array[8]={0,1,2,3,4,5,6,7};
int i,j,k;
int N = 8;
int temp;
j = 0;
for(i = 0; i < N -1; i ++)
{
if(i < j)
{
temp = array[i];
array[i] = array[j];
array[j] = temp;
}
k = N >> 1;
while( k <= j)
{
j = j - k;
k >>= 1;
}
j = j + k;
}
for( i = 0; i < N; i ++)
printf("%d %d\n",i,array[i]);
printf("\n");
return 0;
}
结果:
W r W^r Wr因子的分布
规律:第m级, W 2 m + 1 r W_{2^{m+1}}^r W2m+1r
在fft中,乘法主要来自旋转因子。因为 W r = c o s ( 2 π r N ) − j s i n ( 2 π r N ) , 所 以 在 对 W r W^r=cos(\frac{2{\pi}r}{N})-jsin(\frac{2{\pi}r}{N}),所以在对W^r Wr=cos(N2πr)−jsin(N2πr),所以在对Wr相乘时,必须产生相应的正,余弦函数。编程时,产生正,余弦函数两个方法,一个是在每一步直接产生,二是在程序开始前预先计算W^r存于一个数组中,等效一个正,余弦函数的“表”。
蝶形单元
第m级,有
X m + 1 ( p ) = x m ( p ) + W N r X m ( q ) X_{m+1}(p)=x_{m}(p)+W_N^{r}X_m(q) Xm+1(p)=xm(p)+WNrXm(q)
X m + 1 ( q ) = x m ( p ) − W N r X m ( q ) X_{m+1}(q)=x_{m}(p)-W_N^{r}X_m(q) Xm+1(q)=xm(p)−WNrXm(q)
"组"的概念
每一级的N/2个蝶形单元可以分成若干组,每一组有着相同的结构及W^r因子分布。
第m级的组数是 N / 2 m + 1 N/2^{m+1} N/2m+1.
程序
A.1<<i是把1左移i位,每次左移以为就是乘以2,所以1<<i的结果是1乘以2的i次方 i<<1就是把i左移一位,即i乘以2,假如i=5,最后结果就是5*2=10
B.FFT蝶形算法使用三重循环,下一层的数据都是由上一层计算得到的。
#include <stdio.h>
#include<stdlib.h>
#define pi 3.1415
#define N 64
typedef struct
{
double real,imag;
} complex; //复数结构
complex fft_outa[100];//a的全部数据
complex fft_outao[100];//a的单个数据
complex fft_outb[100];//b的全部数据
complex fft_outbo[100];//b的单个数据
double x_r[N],x_i[N];
double y_r[N],y_[N];
void Rader(double x_r[N],double x_i[N])//雷德算法排序
{
int i,j,k;
j=0;
double t_r,t_i;
for(i=0; i<N-1; i++)
{
if(i<j)
{
t_r=x_r[i];
t_i=x_i[i];
x_r[i]=x_r[j];
x_i[i]=x_i[j];
x_r[j]=t_r;
x_i[j]=t_i;
}
k=N>>1;
while(k<=j)
{
j=j-k;
k>>=1;
}
j=j+k;
}
}
void dit(double x_r[N],double x_i[N])//三层循环
{
int m,r,p,q;
int step,factor_step;
int a,b,i;
double t_real,t_imag;
double w_re,w_im;
for(m=0; m<log2(N); m++) //第一层循环
{
a=1;
for(i=0;i<m+1;i++)
a=a*2;
b=a/2;
w_re=cos(pi/b);
w_im=-sin(pi/b);
step=1<<(m+1);
// factor_step=N>>(m+1);//旋转因素变化速度
//初始化旋转因子
double factor_real=1.0;
double factor_imag=0.0;
for(r=0; r<step/2; r++)
{
// w_re=cos(2*pi*(r+1)*factor_step/N);
// w_im=-sin(2*pi*(r+1)*factor_step/N);
for(p=r; p<N; p+=step)
{
q=p+step/2;//蝶形运算的两个输入
t_real=x_r[q]*factor_real-x_i[q]*factor_imag;
t_imag=x_r[q]*factor_imag+x_i[q]*factor_real;
x_r[q]=x_r[p]-t_real;
x_i[q]=x_i[p]-t_imag;
x_r[p]=x_r[p]+t_real;
x_i[p]=x_i[p]+t_imag;
}
t_real=factor_real*w_re-factor_imag*w_im;
t_imag=factor_real*w_im+factor_imag*w_re;
factor_real=t_real;
factor_imag=t_imag;
}
}
}
int main()
{
int n,i;
double y[N];
for(n=0; n<N; n++)//输入信号
{
x_r[n]=cos(n*pi/6);
x_i[n]=0;
}
Rader(x_r,x_i);
dit(x_r,x_i);
for(i=0; i<N; i++)//输出fft
{
y[i]=sqrt(x_r[i]*x_r[i]+x_i[i]*x_i[i]);
printf("%d %f\n",i,y[i]);
}
}
c结果:
matlab的程序结果:
N=64;
n=0:1:N;
x=cos(n*pi/6);
X=fft(x);
figure
stem(n,abs(X))
对比结果还是有差异。
FFT过程中,对C理解还不够,对蝶形运算算法还不理解,对三重循环还不熟悉 。
修改了,matlab的程序
>> n=0:1:63;
>> x=cos(n*pi/6);
>> y=fft(x);
>> figure
>> stem(n,abs(y))
得出和C一样的结果