【matlab】从DFT到FFT的推导
本文地址https://www.cnblogs.com/jacob1934/p/10478257.html
DFT（离散傅立叶变换）的原理：
要比较两组长度相同的数据的相似性，只需要将两组数据点乘，再求和就行了。
假设两组数据分别为a[N-1:0]和b[N-1:0]，他们的相似性（记为函数r(a,b,N)吧）为
r(a,b,N)=a[0]*b[0] + a[1]*b[1] + a[2]*b[2] + ......+ a[N-2]*b[N-2] + a[N-1]*b[N-1];
如果输入信号为一维数组data，其长度为N，采样率fs=480K。想要求得data的频谱分布，频谱精度为fs/M，也就是输入信号与频率为0*fs/M，1*fs/M，2*fs/M，3*fs/M ...... (M-2)*fs/M，(M-1)*fs/M的波形的相似度。自然是生成长度为N的M组波形，分别与data点乘求和。
那么，测试波形从什么相位开始呢？输入信号中的各个频率也是有相位的，万一相位错了，岂不是都抵消了？！
测试波形干脆用互相正交的cos(x)和sin(x)不就好了！每组抵消部分都能在另一组上体现出来，那怎么算最终结果呢？
cos和sin，那不就是水平轴和垂直轴的么！！
那么，输入信号在该频率上的幅度就是cos和sin的平方根，arctan(sin/cos)就是输入信号的相位了。
等等，有点眼熟，这不就是复数么？
W(k) = cos(k) + sin(k)*i；
一般的，测试波形的相位从0开始，并逐点减少（顺时针旋转），这样就能得到正增长的频谱。而得到的相位和待测波形中sin分量相差90°，或者说与cos分量对齐。
上面说输入信号的各个频率分量都是有相位的，可sin函数值关于0，PI/4，PI/2，PI*3/4都是对称的。我怎么知道，它是相位增加产生的，还是相位减少产生的？
如果输入信号也是由两组互相正交的sin和cos组成的呢？sin(1/8pi)和sin(3/8pi)值相等，cos(1/8pi)和cos(3/8pi)可不相等。
哪来的相位减少的情况？正交下变频中，输入的实数(一维)信号，分别和互相正交(还是cos和sin)的两组本振相乘，得到互相正交的两组信号I和Q。输入信号中频率高于本振的部分被搬移至0~+∞，输入信号的频率低于本振的部分被搬移至-∞~0。这些“负频率”，就是相位减少的。
data也变成复数了，怎么办？还是点乘求和，只不过，是复数和复数的点乘求和。
实数信号能直接和复数测试波形比较么？能啊。实数不正是虚部为0的复数么？不过使用实数信号，频谱将会关于频率0对称，也关于频率fs/2对称。
fs = 480e3; %采样率 N = 2048; %信号长度 M = 2048; %求的DFT点数 %生成待测波形 freq_list=[15e3,-22.5e3,-25e3]; %频率低于本振的信号在正交下变频后，将表现为【负频率】信号 amp_list=[1,0.5,0.5]; src_sin_t=[]; src_cos_t=[]; for i=0:N - 1 src_sin_t(i + 1)=0; src_cos_t(i + 1)=0; for j=1:length(freq_list) freq = freq_list(j); amp = amp_list(j); if(freq > 0) k=i; else k=0-i; end src_sin_t(i + 1) = src_sin_t(i + 1) + amp*sin(2*pi*abs(freq)*k/fs); src_cos_t(i + 1) = src_cos_t(i + 1) + amp*cos(2*pi*abs(freq)*k/fs); end data(i + 1) = complex(src_cos_t(i + 1),src_sin_t(i + 1)); %复数信号 %data(i + 1) = complex(src_cos_t(i + 1),0); %实数信号 end figure(1); subplot(2,1,1); plot(src_sin_t); grid on; subplot(2,1,2); plot(src_cos_t); grid on; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %对测试波形求DFT r_dft=[]; w=[]; %freq=0; phase_inc=0; for i=1:M phase=0; %生成测试波形 for k=1:N radian = 2*pi*phase; %w(i,k) = complex(cos(radian),sin(radian)); w(i,k) = exp(1j*radian); %复指数形式 %phase=phase - freq/fs; phase = phase - phase_inc; end %freq=freq + fs/M; phase_inc = phase_inc + 1/M; %比较待测波形与测试波形 r_dft(i) = 0; for k = 1:N r_dft(i) = r_dft(i) + data(k) * w(i,k); end r_dft(i) = r_dft(i)/N; end %显示 x=[]; q=[]; for i=1:M if(i <= M/2) q(i + M/2) = r_dft(i); else q(i - M/2) = r_dft(i); end %x轴 x(i)=(i - 1 - M/2)*fs/M/1000; end figure(2); subplot(2,1,1); plot(x,abs(q)); grid on; subplot(2,1,2); plot(x,atan2(real(q),imag(q))/pi*180); grid on; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
上面的matlab代码中，M个频率点，每个频率点都要做N个复数乘加运算。
对于频率(k/M*fs)，k是第几个频率点，M是频率点的总数，fs是输入信号的采样频率。测试波形的第i个点的波形值是
freq=k/M*fs; sin_t=sin(2*pi*i*freq/fs);
cos_t=cos(2*pi*i*freq/fs);
约掉fs
sin_t=sin(2*pi*i*k/M);
cos_t=cos(2*pi*i*k/M);
写成复指数形式，exp(-1j*radian)是matlab的语法，就是辐角=radian，模值=1的复数。
w(k,M,i)=exp(-1j*2*pi*i*k/M);
共M个频率点，待测波形长度为N，第k个频率点，是这样运算的
DFT(k,M,N)=
data(0)*exp(-1j*0*k/M*2*pi) +
data(1)*exp(-1j*1*k/M*2*pi)+
data(2)*exp(-1j*2*k/M*2*pi)+
data(3)*exp(-1j*3*k/M*2*pi)+
data(4)*exp(-1j*4*k/M*2*pi)+
data(5)*exp(-1j*5*k/M*2*pi)+
......+
data(N-1)*exp(-1j*(N-1)*k/M*2*pi);
根据【复数相乘，辐角相加，模值相乘】，可以写成这样：
DFT(k,M,N)= data(0)*exp(-1j*0*k/M*2*pi) + data(1)*exp(-1j*0*k/M*2*pi)*exp(-1j*1*k/M*2*pi) +
data(2)*exp(-1j*2*k/M*2*pi) + data(3)*exp(-1j*2*k/M*2*pi)*exp(-1j*1*k/M*2*pi) + data(4)*exp(-1j*4*k/M*2*pi) + data(5)*exp(-1j*4*k/M*2*pi)*exp(-1j*1*k/M*2*pi) + ...... + data(N-2)*exp(-1j*(N-2)*k/M*2*pi) + data(N-1)*exp(-1j*(N-2)*k/M*2*pi)*exp(-1j*1*k/M*2*pi);
为了方便，取M=N=8，提取公共项
DFT(k,8,8)=
(
data(0)*exp(-1j*0*k/8*2*pi) +
data(2)*exp(-1j*2*k/8*2*pi) + // 2/8=1/4，约掉
data(4)*exp(-1j*4*k/8*2*pi) +
data(6)*exp(-1j*6*k/8*2*pi)
)
+
(
data(1)*exp(-1j*0*k/8*2*pi) +
data(3)*exp(-1j*2*k/8*2*pi) +
data(5)*exp(-1j*4*k/8*2*pi) +
data(7)*exp(-1j*6*k/8*2*pi)
)*exp(-1j*1*k/8*2*pi) // 当k>=4时，这里能提一个exp(-1j*1*4/8*2*pi)，也就是旋转半个圆周，等于改变符号
=
(
data(0)*exp(-1j*0*k/4*2*pi) +
data(2)*exp(-1j*1*k/4*2*pi) +
data(4)*exp(-1j*2*k/4*2*pi) +
data(6)*exp(-1j*3*k/4*2*pi)
)
+
(
data(1)*exp(-1j*0*k/4*2*pi)+
data(3)*exp(-1j*1*k/4*2*pi)+
data(5)*exp(-1j*2*k/4*2*pi)+
data(7)*exp(-1j*3*k/4*2*pi)
)*exp(-1j*1*k/8*2*pi)
= // 继续分解
(
(
data(0) + // exp(-1j*0*k/2*2*pi)，等于没乘，直接省略。
data(4)*exp(-1j*1*k/2*2*pi)
)
+
(
data(2) +
data(6)*exp(-1j*1*k/2*2*pi)
)*exp(-1j*1*k/4*2*pi) // k>=4时，这里能提一个exp(-1j*1*4/4*2*pi)，旋转一个圆周，等于没旋转
)
+
(
(
data(1)+
data(5)*exp(-1j*1*k/2*2*pi)
)
+
(
data(3)+
data(7)*exp(-1j*1*k/2*2*pi)
)*exp(-1j*1*k/4*2*pi)
)*exp(-1j*1*k/8*2*pi)
代入k，
DFT(0,8,8) = ( ( data(0) + data(4)*exp(-1j*1*0/2*2*pi) ) + ( data(2) + data(6)*exp(-1j*1*0/2*2*pi) )*exp(-1j*1*0/4*2*pi) ) + ( ( data(1)+ data(5)*exp(-1j*1*0/2*2*pi) ) + ( data(3)+ data(7)*exp(-1j*1*0/2*2*pi) )*exp(-1j*1*0/4*2*pi) )*exp(-1j*1*0/8*2*pi)
DFT(1,8,8) = ( ( data(0) - data(4)*exp(-1j*1*0/2*2*pi) ) + ( data(2) - data(6)*exp(-1j*1*0/2*2*pi) )*exp(-1j*1*1/4*2*pi) ) + ( ( data(1) - data(5)*exp(-1j*1*0/2*2*pi) ) + ( data(3) - data(7)*exp(-1j*1*0/2*2*pi) )*exp(-1j*1*1/4*2*pi) )*exp(-1j*1*1/8*2*pi)
DFT(2,8,8) = ( ( data(0) + data(4)*exp(-1j*1*0/2*2*pi) ) - ( data(2) + data(6)*exp(-1j*1*0/2*2*pi) )*exp(-1j*1*0/4*2*pi) ) + ( ( data(1)+ data(5)*exp(-1j*1*0/2*2*pi) ) - ( data(3)+ data(7)*exp(-1j*1*0/2*2*pi) )*exp(-1j*1*0/4*2*pi) )*exp(-1j*1*2/8*2*pi)
DFT(3,8,8) = ( ( data(0) - data(4)*exp(-1j*1*0/2*2*pi) ) - ( data(2) - data(6)*exp(-1j*1*0/2*2*pi) )*exp(-1j*1*1/4*2*pi) ) + ( ( data(1) - data(5)*exp(-1j*1*0/2*2*pi) ) - ( data(3) - data(7)*exp(-1j*1*0/2*2*pi) )*exp(-1j*1*1/4*2*pi) )*exp(-1j*1*3/8*2*pi)
DFT(4,8,8) = ( ( data(0) + data(4)*exp(-1j*1*0/2*2*pi) ) + ( data(2) + data(6)*exp(-1j*1*0/2*2*pi) )*exp(-1j*1*0/4*2*pi) ) - ( ( data(1) + data(5)*exp(-1j*1*0/2*2*pi) ) + ( data(3) + data(7)*exp(-1j*1*0/2*2*pi) )*exp(-1j*1*0/4*2*pi) )*exp(-1j*1*0/8*2*pi)
DFT(5,8,8) = ( ( data(0) - data(4)*exp(-1j*1*0/2*2*pi) ) + ( data(2) - data(6)*exp(-1j*1*0/2*2*pi) )*exp(-1j*1*1/4*2*pi) ) - ( ( data(1) - data(5)*exp(-1j*1*0/2*2*pi) ) + ( data(3) - data(7)*exp(-1j*1*0/2*2*pi) )*exp(-1j*1*1/4*2*pi) )*exp(-1j*1*1/8*2*pi)
DFT(6,8,8) = ( ( data(0) + data(4)*exp(-1j*1*0/2*2*pi) ) - ( data(2) + data(6)*exp(-1j*1*0/2*2*pi) )*exp(-1j*1*0/4*2*pi) ) - ( ( data(1) + data(5)*exp(-1j*1*0/2*2*pi) ) - ( data(3) + data(7)*exp(-1j*1*0/2*2*pi) )*exp(-1j*1*0/4*2*pi) )*exp(-1j*1*2/8*2*pi)
DFT(7,8,8) = ( ( data(0) - data(4)*exp(-1j*1*0/2*2*pi) ) - ( data(2) - data(6)*exp(-1j*1*0/2*2*pi) )*exp(-1j*1*1/4*2*pi) ) - ( ( data(1) - data(5)*exp(-1j*1*0/2*2*pi) ) - ( data(3) - data(7)*exp(-1j*1*0/2*2*pi) )*exp(-1j*1*1/4*2*pi) )*exp(-1j*1*3/8*2*pi)
整理一下，从下至上的过程
################# ######################################### ################################# #################################
#A(0)= data(0) # #B(0)= A(0) + # #C(0)= B(0) + # #D(0)= C(0) + #
# # # A(1)*exp(-1j*1*0/2*2*pi) # # B(2)*exp(-1j*1*0/4*2*pi)# # C(4)*exp(-1j*1*0/8*2*pi)#
# # # # # # # #
#A(1)= data(4) # #B(1)= A(0) - # #C(1)= B(1) + # #D(1)= C(1) + #
# # # A(1)*exp(-1j*1*0/2*2*pi) # # B(3)*exp(-1j*1*1/4*2*pi)# # C(5)*exp(-1j*1*1/8*2*pi)#
################# ######################################### # # # #
#A(2)= data(2) # #B(2)= A(2) + # #C(2)= B(0) - # #D(2)= C(2) + #
# # # A(3)*exp(-1j*1*0/2*2*pi) # # B(2)*exp(-1j*1*0/4*2*pi)# # C(6)*exp(-1j*1*2/8*2*pi)#
# # # # # # # #
#A(3)= data(6) # #B(3)= A(2) - # #C(3)= B(1) - # #D(3)= C(3) + #
# # # A(3)*exp(-1j*1*0/2*2*pi) # # B(3)*exp(-1j*1*1/4*2*pi)# # C(7)*exp(-1j*1*3/8*2*pi)#
################# ######################################### ################################# # #
#A(4)= data(1) # #B(4)= A(4) + # #C(4)= B(4) + # #D(4)= C(0) - #
# # # A(5)*exp(-1j*1*0/2*2*pi) # # B(6)*exp(-1j*1*0/4*2*pi)# # C(4)*exp(-1j*1*0/8*2*pi)#
# # # # # # # #
#A(5)= data(5) # #B(5)= A(4) - # #C(5)= B(5) + # #D(5)= C(1) - #
# # # A(5)*exp(-1j*1*0/2*2*pi) # # B(7)*exp(-1j*1*1/4*2*pi)# # C(5)*exp(-1j*1*1/8*2*pi)#
################# ######################################### # # # #
#A(6)= data(3) # #B(6)= A(6) + # #C(6)= B(4) - # #D(6)= C(2) - #
# # # A(7)*exp(-1j*1*0/2*2*pi) # # B(6)*exp(-1j*1*0/4*2*pi)# # C(6)*exp(-1j*1*2/8*2*pi)#
# # # # # # # #
#A(7)= data(7) # #B(7)= A(6) - # #C(7)= B(5) - # #D(7)= C(3) - #
# # # A(7)*exp(-1j*1*0/2*2*pi) # # B(7)*exp(-1j*1*1/4*2*pi)# # C(7)*exp(-1j*1*3/8*2*pi)#
################# ######################################### ################################# #################################
所以，计算将分为两步：
【1】按照奇偶重新排列输入数据
【2】从DFT-2开始向上迭代
编写Matlab函数
function [r] = fft_jacob (data,N,invert) %按照奇偶重新排列，从上至下迭代 r=data; for n=log2(N):-1:2 blk_head=0; blk_len=power(2,n); p=[]; for j=1:N/blk_len for k=0:2:blk_len - 2 a=blk_head + k/2; c=blk_head + k; p(a + 1) = r(c + 1); b=blk_head + k/2 + power(2,n - 1); d=blk_head + k + 1; p(b + 1) = r(d + 1); end blk_head = blk_head + blk_len; end r=p; end %从2点DFT开始从下至上迭代 for n=1:log2(N) blk_head = 0; blk_len = power(2,n); p=[]; for j=1:N/blk_len for k=0:blk_len - 1 a=blk_head + rem(k,power(2,n - 1)); b=blk_head + rem(k,power(2,n - 1)) + power(2,n - 1); c=rem(k,power(2,n - 1)); d=power(2,n); if(invert == 1) f=exp(1j*1*c/d*2*pi); else f=exp(-1j*1*c/d*2*pi); end v=blk_head + k; if(k < power(2,n - 1)) p(v + 1) = r(a+1) + r(b+1)*f; else p(v + 1) = r(a+1) - r(b+1)*f; end end blk_head = blk_head + blk_len; end r=p; end if(invert == 0) r=r/N; end
测试用的文件，生成测试波形--计算FFT---对FFT结果计算IFFT---对IFFT结果再求FFT
fs = 480e3; %采样率 N = 2048; %信号长度和FFT点数 M= 2048; %生成待测波形 freq_list=[15e3,-22.5e3,-25e3]; %频率低于本振的信号在正交下变频后，将表现为【负频率】信号 amp_list=[1,0.5,0.5]; src_sin_t=[]; src_cos_t=[]; for i=0:N - 1 src_sin_t(i + 1)=0; src_cos_t(i + 1)=0; for j=1:length(freq_list) freq = freq_list(j); amp = amp_list(j); if(freq > 0) k=i; else k=0-i; end src_sin_t(i + 1) = src_sin_t(i + 1) + amp*sin(2*pi*abs(freq)*k/fs); src_cos_t(i + 1) = src_cos_t(i + 1) + amp*cos(2*pi*abs(freq)*k/fs); end data(i + 1) = complex(src_cos_t(i + 1),src_sin_t(i + 1)); %复数信号 %data(i + 1) = complex(src_cos_t(i + 1),0); %实数信号 end figure(1); subplot(2,1,1); plot(src_sin_t); grid on; subplot(2,1,2); plot(src_cos_t); grid on; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %对测试波形求FFT t0=cputime; r=fft_jacob (data,N,0); t=cputime-t0 %显示 x=[]; q=[]; for i=1:M if(i <= M/2) q(i + M/2) = r(i); else q(i - M/2) = r(i); end %x轴 x(i)=(i - 1 - M/2)*fs/M/1000; end figure(2); subplot(2,1,1); plot(x,abs(q)); grid on; subplot(2,1,2); plot(x,atan2(real(q),imag(q))/pi*180); grid on; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %求IFFT r_ifft=fft_jacob (r,N,1); figure(3); subplot(2,1,1); plot(imag(r_ifft)); grid on; subplot(2,1,2); plot(real(r_ifft)); grid on; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %对r_ifft求FFT q_fft=fft_jacob (r_ifft,N,0); %显示 x=[]; q=[]; for i=1:M if(i <= M/2) q(i + M/2) = q_fft(i); else q(i - M/2) = q_fft(i); end %x轴 x(i)=(i - 1 - M/2)*fs/M/1000; end figure(4); subplot(2,1,1); plot(x,abs(q)); grid on; subplot(2,1,2); plot(x,atan2(real(q),imag(q))/pi*180); grid on;
我写的dft_test.m，2048点，运算耗时10.6S
我写的fft，2048点，运算耗时0.0624S
matlab自带函数fft，2048点，耗时0.0312S
我写的fft函数中有多处可以优化的地方，比如取数组下标都运算了两遍(matlab的数组下标从1开始)，比如p和r的初始化为空，比如power和rem
IFFT可以将频谱和相位转换为波形，与FFT的运算过程相比，测试波形的旋转方向相反。至于IFFT的原理，看看我计算测试波形的过程和IDFT是不是有些相似？
我这里推导的FFT算法，也被称为基2的时间分治FFT算法，是最常见的FFT算法，由J. W. Cooley和 John Tukey 于1965年提出，wikipedia链接Cooley–Tukey FFT algorithm
浙公网安备 33010602011771号