FFT快速傅立叶变换-Matlab和c语言实现

毕业设计做的是基于SOPC的数字频谱分析仪，所以用了一下Matlab。其实数字信号处理学的很差。

FFT来历很复杂，但是他的结果很简单。但凡xxx变换，其实就是一堆数经过运算的到另外一堆数。

对于FFT，输入的数的个数和输出的数的个数相等，这就更妙了。

我们把采样到的数据放入一个虚数数组，数组下标是采样的时间点，计算完成后结果还是存储在这个数组中，

但是下标的含义变成了频率。

Matlab自带的fft函数功能很全很强大，但是这上一个自制的m程序：

function y=myfft(xr,n)

%信号长度为2的整数次方，基于时间抽选的FFT

p=0:n-1; %开始倒位序
nu=log2(n);
p1=p;
b=zeros(1,n);
for t=1:nu;
    p2=floor(p1/2);
    b=b*2+(p1-2*p2);
    p1=p2;
end;
yr(p+1)=xr(b+1);
xr=yr;   %倒位序结束
t=0:n/2-1; %计算因子w开始（只计算w0到w n/2-1）
for v=0:n/2-1;
    w=exp(-2*i*pi*t/n);
end; %计算因子w结束
for m=1:nu;% 计算x（k）开始
    h=2^(m-1);
    k=1;
    while(k<n+1)
        for t=1:h;
            y=bitshift(k-1,nu-m,nu)+1;%求w的幂次数
            xch(k)=xr(k)+w(y)*xr(k+h);
            k=k+1;
        end;
        for t=1:h;
            y=bitshift(k-1-h,nu-m,nu)+1;%求w的幂次数
            xch(k)=xr(k-h)-xr(k)*w(y);
            k=k+1;
        end;
    end;
    xr=xch;
end;%计算x（k）结束
y=xr       %输出变换后的结果

注视就不解释了，现在看function y=myfft(xr,n)

这里xr是虚数，其实就是一个2行n列的矩阵，

更通俗的说，就是2个n长的数组，类型默认是双精度浮点，同样接受返回值的y跟xr一样。

在Matlab下保存这个文件为myfft.m，在同一工作空间就可以用了。

调用这个函数,实现一段信号的幅度谱显示：

%产生两个正弦信号加白噪声
N=2^10;
f1=.1;f2=.2;fs=1;
a1=5;a2=3;
w=2*pi/fs;
x=a1*sin(w*f1*(0:N-1))+a2*sin(w*f2*(0:N-1))+randn(1,N);

subplot(2,2,1);
plot(x(1:N/4)); //只画四分之一长度好看点
title('原始时域信号');

%应用FFT求频谱
x=myfft(x,N);
subplot(2,1,2);
abs_data=abs(x);

for i=1:N
    abs_data(i)=abs_data(i)*2/N;
end
abs_data(1)=abs_data(1)/2;%直流再除以2

f=0:N/2-1; //由于结果的对称性，只显示一半
plot(f,abs_data(1:N/2));
title(‘幅度谱');

%IFFT逆变换，为了验证正变换函数的正确性
y=ifft(x);
subplot(2,2,3);
plot(real(y(1:N/4)));
title('逆变换后时域信号');

由程序就知道了，输入数据为N点的离散时间信号，实部为数据，虚部为0，输出数据也是N点的虚数，

虚数的模反应的就是信号的幅度大小，但是，除了第一个点意外，其余N-1个点的模值都是改频率点幅度值的N/2倍，

第一个点是直流信号幅度的N倍。如果时间信号的采样频率是fs的话，那么幅度谱中的第最高频率是fs/2，也就是第N/2+1点的频率是fs/2，

FFT可以将信号分解得出各个分量的幅度和相位，但是也得考虑时间域的截断造成频谱泄露，使得实际并不是这样，分辨率和分辨力是两个概念。

再上c程序：

void FFT(pJComplex dataIN)

{

uint i,j;

uint mask_1,mask_2;

uint tempIndex_1,tempIndex_2;

uint t,m_1,m_2,n_1;

JComplex buf[LENGTH/2];

 

for(i=0;i<LENGTH;i++)

{

    mask_1 = 1;

    mask_2 = 1<<(LEVEL-1);

    tempIndex_1 = 0;

    for(j=0;j<LEVEL;j++)

    {

      if(i&mask_1)

        tempIndex_1 |= mask_2;

      mask_1 <<= 1;

      mask_2 >>= 1;

    }

    dataIN[i].real = dataIN[tempIndex_1].image;

    dataIN[tempIndex_1].image = 0;

}

for(t=1;t<=LEVEL;t++)

{

    m_1 = 1<<t;

    m_2 = m_1>>1;

    n_1 = LENGTH/m_1;

    for(i=0;i<n_1;i++)

    {

      for(j=0;j<m_2;j++)

      {

        tempIndex_1 = (i<<t)+j;

        tempIndex_2 = tempIndex_1 + m_2;

        Mul(&dataIN[tempIndex_2],&WN[j*n_1],&buf[j]);

        Sub(&dataIN[tempIndex_1],&buf[j],&dataIN[tempIndex_2]);

        Add(&dataIN[tempIndex_1],&buf[j],&dataIN[tempIndex_1]);

      }

    }

}

return;

其中WN是计算好的WN表，Mul等是复数乘法加法减法。。。

板子实际结果与Matlab计算结果：

}

 看图中，我们期望的是只得到3根纯粹的谱线：直流，sin，cos，但是正余弦分量呈现为塔状，如果正弦和余弦分量的频率值比较接近，那么彼此交叠的现象就严重了，

这就是频谱泄漏，怎么改善呢？可以延长计算点数，使其更加接近理论上的无限长，但是计算点数的延长使得算法耗时增加。

实物：板子芯片是EP2C8Q208C8