FFT原理及实现

目录

1.基础介绍

2.FT介绍

3.DFT计算

4.FFT计算

5.DFT和FFT仿真

6.完整matlab工程代码

---

 1.基础介绍

由傅里叶变换可知：

1. 时域离散---->频域周期
2. 时域连续---->频域非周期
3. 时域周期---->频域离散
4. 时域非周期---->频域连续

一般通信系统信号经过AD转换后都为离散非周期信号，因此对应的频域则为连续周期信号。

连续信号无法使用计算机或其他数字处理器进行处理，所以连续频域信号也需要离散化，即对频域进行抽样，频域信号与抽样频域函数相乘，所对应的时域运算就是和抽样序列进行卷积。【时域相乘---->频域卷积，频域相乘---->时域卷积】

2.FT介绍

由上图可知：

时域:模拟信号.*采样序列---->采样后信号【ADC模数变换】【.*=乘，*=卷积】

图1.*图3=图5

频域:模拟信号频谱*采样序列频谱---->采样后频谱【DTFT离散时间傅里叶变换】

图2*图4=图6

频域:采样后频谱.*频域采样序列频谱---->离散采样后频谱【DFT离散傅里叶变换】

图6.*图8=图10

时域:采样后信号*频域采样序列【IDFT离散傅里叶逆变换】

图5*图7=图9；图10---->IDFT---->图9。

针对DFT运算进行计算效率优化和复杂度优化，提出FFT运算，其本质仍为DFT。

3.DFT计算

傅里叶变换对：

对时间t->n和频域f->k上使用有限数量∞->N的采样[积分->求和]即DFT得到：

IDFT为：

DFT运算举例：[乘加运算量/复杂度=O(N^2)]

设有N=3个采样点数据; ; ;

根据DFT公式有：

k=0：

k=1：

k=2：

4.FFT计算

注：FFT方便处理长度N=2^P的情况，如果长度不是2的整数次幂，则补零即可。

首先将长度为N的序列x(n)[n=0…N-1]分为两个奇偶序列之和：x1(m)=x(2n); x2(m)=x(2n+1)；x1和x2两个序列的的长度都为M=N/2。

根据DFT公式有：

分成奇偶两部分有：[k=0…N-1]

又有[k=0…M-1]

因此：

 

因此：

以此类推可以将x1(m)和x2(m)不断的进行奇偶分解，

[乘加运算量/复杂度=O(Nlog2^N)]，例如8个采样点的数据：

0,1,2,3,4,5,6,7-->0,2,4,6;1,3,5,7-->0,4;2,6;1,5;3,7-->0;4;2;6;1;5;3;7。0级-->1-->2-->3级。

将分解过程分步定义，建立递归过程X[k=0-7]-->X`[k=0-3]-->X``[k=0-1]-->X```[k=0]；

3级有：---->  [最后第3级级k固定为0]

2级有：---->

1级有：

0级有：

5.DFT和FFT仿真

1、设计一个以1M采样率的100KHz频率的正弦信号sin_f0和一个200KHz的正弦信号sin_f1，将两者相加得到包含两个频率的信号y，采样点数1024个点。

2、使用matlab自带的fft函数对y进行FFT运算，观察其频谱，在100KHz处和200KHz处各有一个尖峰。

3、实现一个DFT程序，运算长度为128即分辨率=1M/128=7.8KHz，设计完成后对信号y进行运算，输出结果观察其频谱。

4、实现一个FFT程序，运算长度为128即分辨率=1M/128=7.8KHz，设计完成后对信号y进行运算，输出结果观察其频谱，与使用matlab自带函数进行对比。

6.完整matlab工程代码

%------------------------------------
%------------------------------------
%------------------------------------
close all;
clc;
%------------------------------------
%------------------------------------
%------------------------------------
sin_f0=0.1E6;%构造一个100KHz频率的正弦信号y0
sin_f1=0.2E6;%构造一个200KHz频率的正弦信号y1
Fs=1E6;%采样率为1MSPS
N=2^10;
n=0:N-1;
time=n/Fs;
y0=sin(2*pi*sin_f0*time);
y1=sin(2*pi*sin_f1*time);
y=y0+y1;%将100KHz的信号200KHz的信号相加，生成包含两个频率分量的原始信号
%------------------------------------
%------------------------------------
%------------------------------------
base_len = length(y);
f_lab0 = zeros(1,base_len);%将采样点换算频域横坐标轴--频率
for i = 1:base_len
    f_lab0(i) = i*Fs/base_len;
end
f_lab1 = zeros(1,128);%将采样点换算频域横坐标轴--频率
for i = 1:128
    f_lab1(i) = i*Fs/128;
end
%------------------------------------
%------------------------------------
%------------------------------------
base_sig = y;
base_fft = abs(fft(base_sig(1:128)));
subplot(4,1,1);
plot(base_sig(1:1000),'b');grid;%观察原始信号的时域波形
title('原始时域信号');
subplot(4,1,2);
plot(f_lab1,base_fft,'b');grid;%观察原始信号的频域波形
title('工具自带FFT函数计算的频域信号');
%------------------------------------
%------------------------------------
%------------------------------------
dft_in = zeros(1,128);
dft_out = zeros(1,128);
for i = 1:128
    dft_in(i) = base_sig(i);
end

for j=1:128
    dft_out(j) = 0;
    for i=1:128
        dft_out(j) = dft_out(j) + dft_in(i)*complex(cos(2*pi*i*j/128),-sin(2*pi*i*j/128));
    end 
end 

subplot(4,1,3);
plot(f_lab1,abs(dft_out),'b');grid;%观察自实现DFT运算后的频域波形
title('自实现DFT计算的频域信号');
%------------------------------------
%------------------------------------
%------------------------------------
fft_in0 = zeros(1,128); %进行奇偶拆分
fft_in1 = zeros(1,128); 
fft_in2 = zeros(1,128); 
fft_in3 = zeros(1,128); 
fft_in4 = zeros(1,128); 
fft_in5 = zeros(1,128); 
fft_in6 = zeros(1,128);
for i = 1:128
    fft_in0(i) = base_sig(i);
end

for i = 1:64
    fft_in1(i) = fft_in0(2*i-1);
    fft_in1(64+i) = fft_in0(2*i);
end 

for j = 1:2 
    k=(j-1)*64;
    for i = 1:32
        fft_in2(i+k) = fft_in1(2*i-1+k);
        fft_in2(32+i+k) = fft_in1(2*i+k);
    end 
end 

for j = 1:4 
    k=(j-1)*32;
    for i = 1:16
        fft_in3(i+k) = fft_in2(2*i-1+k);
        fft_in3(16+i+k) = fft_in2(2*i+k);
    end 
end 

for j = 1:8 
    k=(j-1)*16;
    for i = 1:8
        fft_in4(i+k) = fft_in3(2*i-1+k);
        fft_in4(8+i+k) = fft_in3(2*i+k);
    end 
end

for j = 1:16 
    k=(j-1)*8;
    for i = 1:4
        fft_in5(i+k) = fft_in4(2*i-1+k);
        fft_in5(4+i+k) = fft_in4(2*i+k);
    end 
end

for j = 1:32 
    k=(j-1)*4;
    for i = 1:2
        fft_in6(i+k) = fft_in5(2*i-1+k);
        fft_in6(2+i+k) = fft_in5(2*i+k);
    end 
end
%------------------------------------
%------------------------------------
%------------------------------------
fft_out0 = zeros(1,128); %进行FFT运算
fft_out1 = zeros(1,128); 
fft_out2 = zeros(1,128); 
fft_out3 = zeros(1,128); 
fft_out4 = zeros(1,128); 
fft_out5 = zeros(1,128); 
fft_out6 = zeros(1,128); 
for i = 1:64
    fft_out0(2*i-1:2*i) = fft_sub_fun(fft_in6(2*i-1),fft_in6(2*i),1,2);
end

for i = 1:32
    fft_out1(4*i-3:4*i) = fft_sub_fun(fft_out0(4*i-3:4*i-2),fft_out0(4*i-1:4*i),2,4);
end 

for i = 1:16
    fft_out2(8*i-7:8*i) = fft_sub_fun(fft_out1(8*i-7:8*i-4),fft_out1(8*i-3:8*i),4,8);
end 

for i = 1:8
    fft_out3(16*i-15:16*i) = fft_sub_fun(fft_out2(16*i-15:16*i-8),fft_out2(16*i-7:16*i),8,16);
end 

for i = 1:4
    fft_out4(32*i-31:32*i) = fft_sub_fun(fft_out3(32*i-31:32*i-16),fft_out3(32*i-15:32*i),16,32);
end 

for i = 1:2
    fft_out5(64*i-63:64*i) = fft_sub_fun(fft_out4(64*i-63:64*i-32),fft_out4(64*i-31:64*i),32,64);
end 

for i = 1:1
    fft_out6 = fft_sub_fun(fft_out5(1:64),fft_out5(65:128),64,128);
end 

subplot(4,1,4);
plot(f_lab1,abs(fft_out6),'b');grid;%观察自实现FFT运算后的频域波形
title('自实现FFT计算的频域信号');

//fft_sub_fun函数代码//

%------------------------------------
%------------------------------------
%------------------------------------
function fft_out = fft_sub_fun(fft_in0,fft_in1,fft_len,N)

fft_out = zeros(1,fft_len*2); 
for i=1:fft_len 
    fft_out(i) = fft_in0(i)+fft_in1(i)*complex(cos(2*pi*i/N),-sin(2*pi*i/N));
    fft_out(i+fft_len) = fft_in0(i)-fft_in1(i)*complex(cos(2*pi*i/N),-sin(2*pi*i/N));
end