数字信号处理知识点总结（三）：离散傅里叶变换（DFT）
本篇文章主要介绍离散傅里叶变换（DFT）的概念、推导、实现及性质以及用DFT实现的线性卷积。
目录
-
-
1. 离散傅里叶变换（DFT）概念
-
2. 离散傅里叶变换（DFT）的公式推导
-
3. 离散傅里叶变换（DFT）的实现
-
4. 离散傅里叶变换（DFT）的性质
-
- 4.1循环反转
- 4.2对称性
- 4.3循环移位
- 4.4循环卷积
-
5. 离散傅里叶变换（DFT）实现的线性卷积
-
1. 离散傅里叶变换（DFT）概念
根据数字信号处理知识点总结（二）中的公式（5）可知，一个非周期离散无限长信号进行离散时间傅里叶变换（DTFT）后频谱是连续周期的，以相距\omega = 2\pi /N的等分点对该频谱进行采样就会产生一个DFS序列（其是离散周期的，如数字信号处理知识点总结（二）的公式（7）所示），对其进行IDFS（离散傅里叶合成）则会产生周期离散的序列\tilde x(n)。
但是，实际上大多数用于计算机处理的信号都不是周期的，而且计算机也只能处理有限长的信号。理论上，有两种方法可以实现：第一种就是在该有限长序列的两边添上无限的零信号，便能运用DTFT的相关理论进行求解；第二种便是对有限长序列进行周期延拓，运用离散时间傅里叶分解（DFS）进行求解。两种方法都是把有限长信号转变为无限长信号进行求解。第一种方法由于需要存储无限多（或者趋向无限多）的零，对计算机存储性能较大，故不允采用。第二种方法中，周期延拓后的周期信号的主值区间（周期序列延拓的第一个周期）就是原本的有限长信号，然后对这个周期信号应用DFS，DFS后得到结果的主值周期 可以定义为离散傅里叶变换（DFT） 。
2. 离散傅里叶变换（DFT）的公式推导
定义一个有限长序列x(n)，它在0 \le n \le N - 1上有N个样本，作为一个N点序列。令\tilde x(n)是用这个N点序列x(n)创建的一个周期为N的周期信号，如公式（1）所示：
\tilde x(n) = \sum\limits_{r = - \infty }^\infty {x(n - rN)} \tag{1}
该式子通过模N运算（modulo-N）可以简化得到公式（2）：
\tilde x(n) = x(n\bmod N)\tag{2}
模N运算可以简单地说明：如果宗量n是在0和N-1之间，那就是它自己；否则，从n开始加或减N的倍数，直到结果是在0和N-1之间为止。x(n)的长度应当是N或小于N时才成立，否则若n取值大于N，则会发生频谱的混叠。公式（2）可以使用公式（3）来表示：
x{((n))_n}\mathop = \limits^\Delta x(n\bmod N)\tag{3}
综上x(n)和\tilde x(n)之间的关系是
\begin{array}{l} \tilde x(n) = x{((n))_N}\\ \tag{4} {\rm{x}}(n) = \tilde x(n){R_N}(n) \end{array}
其中，公式（4）的第一个式子是周期延拓，第二个式子是时窗运算，{R_N}(n)是时窗函数。
由{\tilde X}(k)\mathop = \limits^\Delta DFS[{\tilde x}(n)],可以推导得到一个N点序列的离散傅里叶变换为：
{\rm{X}}(k)\mathop = \limits^\Delta DFT[x(n)] = \left\{ \begin{array}{l} {\rm{\tilde X}}(k),0 \le k \le N - 1\\ 0,other \end{array} \right. = {\rm{\tilde X}}(k){R_N}(k)\tag{5}
同理，由{\tilde x}(n)\mathop = \limits^\Delta IDFS[{\tilde X}(k)]可以推导得到一个N点的DFTX(k)的逆离散傅里叶变换为：
{\rm{x}}(n)\mathop = \limits^\Delta IDFT[X(k)] = \tilde x(n){R_N}(n) \tag{6}
整理公式（5）和公式（6）得到DFT公式：
X(k) = \sum\limits_{n = 0}^{N - 1} {x(n)\exp ( - j\frac{{2\pi }}{N}nk),0 \le k \le N - 1}\tag{7}
以及IDFT公式：
x(n) = \frac{1}{N}\sum\limits_{k = 0}^{N - 1} {X(k)\exp (j\frac{{2\pi }}{N}nk),0 \le n \le N - 1}\tag{8}
3. 离散傅里叶变换（DFT）的实现
以下是使用矩阵相乘的结构对DFT的matlab实现代码：
function [ Xk ] = dft(xn,N)
%complete discrete fourier series coefficents
% %
% [Xk] = dft(xn,N)
% Xk = DFS coeff.array over <= k <= N-1
% xn = One period of periodic signal over 0 <= n <= N-1
% N = Fundamental period of xn
n = [0:1:N-1];
k = [0:1:N-1];
WN = exp(-j*2*pi/N);
nk = n' * k;
WNnk = WN .^nk;
Xk = xn(1:1:N) * WNnk;
end以及相应的IDFT实现，实质上是在以上代码的结构上进行了非常少的改变：
function [ xn] = idft(Xk,N)
%complete inverse discrete fourier series
% %
% [ xn] = idfs(Xk,N)
% Xk = DFS coeff.array over <= k <= N-1
% xn = One period of periodic signal over 0 <= n <= N-1
% N = Fundamental period of xn
n = [0:1:N-1];
k = [0:1:N-1];
WN = exp(-j*2*pi/N);
nk = n' * k;
WNnk = WN .^(-nk);
xn = (Xk(1:1:N) * WNnk)/N;
end4. 离散傅里叶变换（DFT）的性质
本节主要介绍几个DFT特有且重要的性质，包括循环反转、对称性、循环移位、循环卷积。其中，循环反转、循环移位、循环卷积都是因为DFT定义的范围限定在0 \le n \le N - 1而衍生出来的性质。
4.1循环反转
定义循环反转为：
x{(( - n))_N} = \left\{ \begin{array}{l} x(0),n = 0\\ x(N - n),1 \le n \le N - 1 \end{array} \right.\tag{8}
为了将它形象化，想象将序列x(n) 以逆时针方向放在一个圆上， 并使n =0和n = N 重合，那么x( ( - n) )_N 就能看作是以顺时针方向将x(n)放在圆上，所以也称为圆周反转 。在MATLAB 中可用x= x(mod( - n ,N) + 1 ）实现循环反转。在matlab中宗量n从一开始，故要加上mod()函数后要加上1。循环反转的DFT 给出为：
DFT[x{(( - n))_N}] = X{(( - k))_N} = \left\{ \begin{array}{l} X(0),k = 0\\ X(N - k),1 \le k \le N - 1 \end{array} \right..\tag{9}
在matlab中某个直观的例子如下：
4.2对称性
共轭对称：
DFT[{x^*}(n)] = {X^*}{(( - k))_N}\tag{10}
若 x(n)是实序列，则有x(n)={x^*}(n)，利用公式（10）有以下的循环（或圆周）共轭对称性质：
X(k) = {X^*}{(( - k))_N}\tag{11}
进而可以得到：
对于一组离散的N点实序列x(n)，可以分解为循环偶分量 和循环奇分量 ：
\begin{array}{l} {x_{ec}}\mathop = \limits^\Delta \frac{1}{2}[x(n) + x{(( - n))_N}] = \left\{ \begin{array}{l} x(0),n = 0\\ \frac{1}{2}[x(n) + x(N - n)],1 \le n \le N - 1 \end{array} \right.\\ {x_{oc}}\mathop = \limits^\Delta \frac{1}{2}[x(n) - x{(( - n))_N}] = \left\{ \begin{array}{l} 0,n = 0\\ \frac{1}{2}[x(n) - x(N - n)],1 \le n \le N - 1 \end{array} \right. \end{array}\tag{12}
且有：
\begin{array}{l} DFT[{x_{ec}}(n)] = {\mathop{\rm Re}\nolimits} [X(k)] = {\mathop{\rm Re}\nolimits} [X{(( - k))_N}]\\ DFT[{x_{oc}}(n)] = {\mathop{\rm Im}\nolimits} [X(k)] = {\mathop{\rm Im}\nolimits} [X{(( - k))_N}] \end{array}\tag{13}
实现代码：
function [xec,xoc] = circevod(x)
%signal decomposition into circular-even and circular - odd parts
% -----------------------------------------------------
if any(imag(x) ~= 0)
error ('x i s not a real sequence')
end
N = length(x); n = 0: (N - 1);
xec = 0.5*(x+x(mod(-n,N)+1));xoc = -0.5*(x-x(mod(-n,N)+1));
end4.3循环移位
对于循环移位，首先得把x(n)延拓为\tilde x(n)，然后再移m个样本，得到
\tilde x(n - m) = x{((n - {\rm{m}}))_N}\tag{14}
这称为\tilde x(n)的周期移位 。然后再将周期移位转换为N 点序列。所得序列:
\tilde x(n - m){R_N}(n) = x{((n - {\rm{m}}))_N}{R_N}(n)\tag{15}
称为x(n)的循环移位 。可以设想将序列x(n)放在一个圆上， 现在将这个圆旋转m个样本， 并从0 \le n \le N - 1内展开这个序列。它的DFT 给出为：
\begin{array}{l} DFT[x{((n - {\rm{m}}))_N}{R_N}(n)] = W_N^{km}X(k)\\ W_N^{km} = \exp ( - j\frac{ {2\pi }}{N}km) \end{array}\tag{16}
实现代码：
function [y] = cirshftt(x,m,N)
%实现循环移位
% -----------------------------------------------------
% [y] = cirshftt(x,m,N)
% y = output sequence containing the circular shift
% x = input sequence of length <= N
% m = sample shift
% N = size of circular buffer
% Method, y(n) = x( (n-m) mod N)
%Check for length of x
if length(x) > N
error('N must be >= the length of x')
end
x = [x zeros(1,N-length(x))];
n = [0:1:N-1];n = mod(n-m,N);y=x(n+1);
end循环移位例子：
4.4循环卷积
时域方法循环卷积实现代码：
function [y] = circonvt(x1,x2,N)
%N-point c ircular convolution between xl and x2: ( time - domain)
%－－－－－－－－－－－－－－－－－－－－－－－－一－一－－－－－－－－－－－－－－－－－－
% [ y] = circonvt (xl ,x2 ,N)
% y = output sequence containing the circular convol ution
% xl = input sequence of length Nl < = N
% x2 = input sequence of length N2 < = N
% N = size of circular buffer
% Method: y(n) = sum (xl(m)*x2((n - m) mod N)
% Check for length of xl
if length(x1) > N
error('N must be > = the length of xl')
end
% Check for length of x2
if length(x2) > N
error('N must be > = the length of x2')
end
x1 = [x1 zeros(1,N- length( x1))];
x2 = [x2 zeros(1, N - length( x2))];
n = [0:1:N-1];x2 = x2(mod(-n,N)+1);H=zeros(N,N);
for n=1:1:N
H(n,:)=cirshftt(x2,n-1,N);%n从1开始数
end
y = x1*conj(H');
end注意代码的倒数第二行，要乘卷积核H的共轭 才能得到正确的结果。频域方法实现的循环卷积采用IDFT[DFT[x1]*DFT[x2]]的方法即可得到。
5. 离散傅里叶变换（DFT）实现的线性卷积
通过对比循环卷积与线性卷积（详细推导可见《数字信号处理matlab版》中p151-153），当取N=max({N_1},{N_2})作循环卷积时，那么前(M-1)个样本在误差中（也就是跟线性卷积的不一样），这里M=min({N_1},{N_2})。由此对数据进行分块的批处理方法很有用处。这里主要介绍批处理中的重叠保留法 来实现线性卷积：
算法步骤描述
时域方法实现的matlab代码：
function [ y ] = ovrlpsav(x,h,N)
%Overlap-Save method of block convolution
% 此处显示详细说明
% [ y ] = ovrlpsav(x,h,N)
% y = output sequence
% x = input sequence
% h = impul se response
% N = block length
Lenx = length(x);M = length(h);M1 = M - 1;L = N-M1;
h = [h zeros(1,N-M)]; %长度为N
%
x = [zeros(1,M1),x,zeros(1,N-1)]; %zeros(1,N-1)分段最后一段最少只有一个x的样本值
K = floor((Lenx+M1-1)/L);
Y = zeros(K+1,N);
% convolution wi th succesive blocks
for k=0:K
xk = x(k*L+1:k*L+N); %取每一小段x，长N
Y(k+1,:) = circonvt(xk,h,N);
end
Y = Y(:,M:N)';
y = (Y(:))';
end以上是时域方法的线性卷积实现，下面使用频域方法DFT实现线性卷积：
function [y] = dft_ovrlpsav( x,h,N)
%High-speed O飞rerlap-Save method of block convolutions using FFT
% －－－－－－－一－ －－－－－－－－－－ －一－ －－ － －－－－－－－－－－－－－－－
%［y]＝ hsolpsav(x,h,N)
% y = output sequence
% x = input sequence
% h = impulse r esponse
% N = block length (must be a power of two)
N = 2^(ceil(log10(N)/log10(2)));
Lenx = length(x);M = length(h);M1 = M - 1;L = N-M1;
h = dft(h,N); %先做好FFT变换
%
x = [zeros(1,M1),x,zeros(1,N-1)];
K = floor((Lenx+M1-1)/L);
Y = zeros(K+1,N);
% convolution wi th succesive blocks
for k=0:K
xk = dft(x(k*L+1:k*L+N));
Y(k+1,:) = idft(xk.*h);
end
Y = Y(:,M:N)';
y = (Y(:))';
end其中dft()与idft()函数的实现见本文的第三小节。
1.[美]纳•K•英格尔，[美]约翰•G•普罗克斯. 数字信号处理（MATLAB版）（第3版）[M].西安交通大学出版社，2013.7
2.Oppenheim, A. V., Willsky, A. S.,Nawab, S. H. 信号与系统（第二版）[M]. 电子工业出版社， 2013.1
3.John G. Proakis, Dimitris G. Manolakis. 数字信号处理——原理、算法与应用（第四版）[M]. 电子工业出版社，2014.8
4.Steven W. Smith. 实用数字信号处理[M]. 人民邮电出版社，2010.12