快速傅里叶变换FFT的代码实现

最新推荐文章于 2024-05-24 11:22:55 发布

古楼望月

最新推荐文章于 2024-05-24 11:22:55 发布

阅读量1.5k

点赞数 4

分类专栏：数字信号处理文章标签：快速傅里叶变换算法

本文链接：https://blog.csdn.net/YDY5659150/article/details/103436889

版权

数字信号处理专栏收录该内容

9 篇文章 2 订阅

订阅专栏

本文将根据离散傅里叶变换的快速算法编写代码

1，首先列出DFT的变换式和快速计算公式

$X(k)=\sum_{n=0}^{N-1}x(n)W_N^{nk}$

$X(k)=X_1(k)+X_2(k)W_N^k$

$X(k)=X_1(k)-X_2(k)W_N^k$

2，假定N=2^L，那么采用基拆分具有如下树状规律

3，下面在表格里面列出快速算法的一般规律：

	左子树	右子树	序列分组	组内分对	W(N,k)
第1层	b[L-1]=0	b[L-1]=1	2^(L-1)	2^(1-1)	N=2^1,k=[0~2^(1-1)-1]
第2层	b[L-2]=0	b[L-2]=1	2^(L-2)	2^(2-1)	N=2^2,k=[0~2^(2-1)-1]
...	...	...	...	...	...
第i层	b[L-i]=0	b[L-i]=1	2^(L-i)	2^(i-1)	N=2^i,k=[0~2^(i-1)-1]
...	...	...	...	...
第L层	b[L-L]=0	b[L-L]=1	2^(L-L)	2^(L-1)	N=2^L,k=[0~2^(L-1)-1]

4，运算的时候，从序号的高位进行，向低位滑行。需要注意的是取k值，即是蓝色区域的数值需要按位逆序运算。完整的代码如下：

function X_k=fft2_dit(x_n)
N1 =length(x_n);
L =ceil(log2(N1));
N =2^L;
X_t=zeros(1,N);
%如果序列不是2整数次幂，那么进行补零处理
for m=1:N
if(m<=N1)
X_t(m)=x_n(m);
else
X_t(m)=0;
end
end
for m=1:L %层级循环
N =2^m;%W(N,k)中N值
k_max =N/2 -1;%W(N,k)中最大k值
g_max =2^(L-m);%一层有多少组独立的蝶形图
for g=1:g_max %层内分组
for p=0:k_max %组内分对
k=0; %k值逆序运算
for n=1:m-1
k=bitshift(k,1) + bitget(p,n);
end
num_1 =(2*p+0)*g_max+ g;
num_2 =(2*p+1)*g_max+ g;
X1 =X_t(num_1) + X_t(num_2)*exp(-1j*2*pi*k/N);
X2 =X_t(num_1) - X_t(num_2)*exp(-1j*2*pi*k/N);
X_t(num_1)=X1;
X_t(num_2)=X2;
end
end
end

X_k=zeros(1,N);
for m=0:2^L-1
k=0; %输出逆序运算
for n=1:L
k=bitshift(k,1) + bitget(m,n);
end
X_k(m+1)=X_t(k+1);
end
end

5，下面测试按时间基2拆分的正确性

x_n =[0 1 2 3 4 5 6 7];
Xk =fft2_dit(x_n);%测试输出
Xp =fft(x_n);%使用matlab自带的函数
Xk =Xk.'
Xp =Xp.'

测试结果如下
Xk =

28.0000
-4.0000 + 9.6569i
-4.0000 + 4.0000i
-4.0000 + 1.6569i
-4.0000
-4.0000 - 1.6569i
-4.0000 - 4.0000i
-4.0000 - 9.6569i

Xp =

28.0000
-4.0000 + 9.6569i
-4.0000 + 4.0000i
-4.0000 + 1.6569i
-4.0000
-4.0000 - 1.6569i
-4.0000 - 4.0000i
-4.0000 - 9.6569i

6，首先列出DFT的变换式和快速计算公式

$X(k)=\sum_{n=0}^{N-1}x(n)W_N^{nk}$

$x_1(n)=x(n)+x(n+\frac{N}{2})$

$x_1(n)=[x(n)-x(n+\frac{N}{2})]W_N^n$

7,分析之后可以发现，每进行一层蝶形计算之后，序列对折一半。显然对折的比特位从最高位开始，直到最低位结束，层层依次迭代。因为每进行一层运算，频域值得序号进行一次奇偶交换，等于频域值对折的比特位从最低位开始，直到最高位结束。需要注意的是k值不会有逆序的操作。最后得到的频域值，进行一次逆序操作。程序如下：

function X_k=fft2_dif(x_n)
N1 =length(x_n);
L =ceil(log2(N1));
N =2^L;
X_t=zeros(1,N);
%如果序列不是2整数次幂，那么进行补零处理
for m=1:N
if(m<=N1)
X_t(m)=x_n(m);
else
X_t(m)=0;
end
end
for m=1:L %层级循环
N =2^(L-m+1);%W(N,k)中N值
k_max =N/2-1;%W(N,k)中最大k值
g_max =2^(m-1);%一层有多少组独立的蝶形图
for g=0:g_max-1 %层内分组
for k=0:k_max %组内分对
num_1 =(2*g+0)*2^(L-m)+ k+1;
num_2 =(2*g+1)*2^(L-m)+ k+1;
X1 =X_t(num_1) + X_t(num_2);
X2 =(X_t(num_1) - X_t(num_2))*exp(-1j*2*pi*k/N);
X_t(num_1)=X1;
X_t(num_2)=X2;
end
end
end

X_k=zeros(1,N);
for m=0:2^L-1
k=0; %输出逆序运算
for n=1:L
k=bitshift(k,1) + bitget(m,n);
end
X_k(m+1)=X_t(k+1);
end
end

8，下面测试按频率基2拆分的正确性

x_n =[0 1 2 3 4 5 6 7 8 7 6 5 4 3 2 1];
Xk =fft2_dit(x_n);%测试输出
Xp =fft(x_n);%使用matlab自带的函数
Xk =Xk.'
Xp =Xp.'

测试结果如下

Xk =

64.0000
-26.2741 - 0.0000i
0
-3.2398 - 0.0000i
0
-1.4465 + 0.0000i
0
-1.0396 + 0.0000i
0
-1.0396 + 0.0000i
0
-1.4465 + 0.0000i
0
-3.2398 + 0.0000i
0
-26.2741 - 0.0000i

Xp =

64.0000
-26.2741
0
-3.2398
0
-1.4465
0
-1.0396
0
-1.0396
0
-1.4465
0
-3.2398
0
-26.2741

自此，快速傅里叶算法研究到此结束，如果对你有所帮助，请点一个赞，谢谢。

古楼望月

关注

4
点赞
踩
5

收藏

觉得还不错? 一键收藏
0
评论
快速傅里叶变换FFT的代码实现

本文将根据离散傅里叶变换的快速算法编写代码1，首先列出DFT的变换式和快速计算公式2，假定N=2^L，那么采用基拆分具有如下树状规律3，下面在表格里面列出快速算法的一般规律：左子树右子树序列分组组内分对 W(N,k) 第1层 b[L-1]=0 b[L-1]=1 2^(L-1) 2^(1-1) ...
复制链接

扫一扫