DFT

公式：

Ak=∑n=0N−1e−i2πNknan,k∈{0,1,...N−1}A_k = \sum_{n=0}^{N-1}e^{-i\frac{2\pi}{N}kn}a_n , k\in \{0,1,... N-1\} Ak=n=0∑N−1e−iN2πknan,k∈{0,1,...N−1}

一般我们通常写成
Ak=∑n=0N−1WNknan,k∈{0,1,...N−1}A_k = \sum_{n=0}^{N-1}W_{N}^{kn}a_n , k\in \{0,1,... N-1\} Ak=n=0∑N−1WNknan,k∈{0,1,...N−1}

其中

WN=e−i2πNW_{N} = e^{-i\frac{2\pi}{N}} WN=e−iN2π

WNkW_{N}^{k}WNk称为旋转因子, twiddle factor，之所以称之为旋转因子，是因为它在单位圆上根据k值不同而做不同角度的旋转。 (更多阅读请参考欧拉公式和傅里叶变换原理)

下面是N=2,N=4,N=8的旋转，

从图中可以看到，旋转因子具有如下性质

WN0=1WNN2=−1WNk=WNk+NW_{N}^{0}=1 \\ W_{N}^{\frac{N}{2}}=-1 \\ W_{N}^{k}=W_{N}^{k+N} WN0=1WN2N=−1WNk=WNk+N

iDFT

公式：
an=1N∑k=0N−1WN−knAk,k∈{0,1,...N−1}a_n = \frac{1}{N}\sum_{k=0}^{N-1}W_N^{-kn}A_k, k\in \{0,1,... N-1\} an=N1k=0∑N−1WN−knAk,k∈{0,1,...N−1}

具体示例

2个sample的DFT

旋转因子

W2=e−i2π/2=e−iπ=cosπ−isinπ=−1W_2 = e^{-i2\pi/2} = e^{-i\pi} = cos\pi-isin\pi = -1 W2=e−i2π/2=e−iπ=cosπ−isinπ=−1

AkA_kAk:
Ak=∑n=01(−1)knan=(−1)0ka0+(−1)1ka1A_k = \sum_{n=0}^1 (-1)^{kn}a_n = (-1)^{0k}a_0+(-1)^{1k}a_1 Ak=n=0∑1(−1)knan=(−1)0ka0+(−1)1ka1

所以
A0=a0+a1A1=a0−a1A_0 = a_0+a_1 \newline A_1 = a_0-a_1 A0=a0+a1A1=a0−a1

4个sample的DFT

旋转因子

W4=e−i2π/4=e−iπ/2=−iW_4 = e^{-i2\pi/4} = e^{-i\pi/2} = -i W4=e−i2π/4=e−iπ/2=−i

AkA_kAk:
Ak=∑n=03(−1)knan=(−i)0ka0+(−i)1ka1+(−i)2ka2+(−i)3ka3A_k = \sum_{n=0}^3 (-1)^{kn}a_n = (-i)^{0k}a_0+(-i)^{1k}a_1+(-i)^{2k}a_2+(-i)^{3k}a_3 Ak=n=0∑3(−1)knan=(−i)0ka0+(−i)1ka1+(−i)2ka2+(−i)3ka3

所以
A0=a0+a1+a2+a3A1=a0−ia1−a2+ia3A2=a0−a1+a2−a3A3=a0+ia1−a2−ia3A_0 = a_0+a_1+a_2+a_3 \newline A_1 = a_0-ia_1-a_2+ia_3 \newline A_2 = a_0-a_1+a_2-a_3 \newline A_3 = a_0+ia_1-a_2-ia_3 A0=a0+a1+a2+a3A1=a0−ia1−a2+ia3A2=a0−a1+a2−a3A3=a0+ia1−a2−ia3

可以把它抽取出来参数，作为一个矩阵相乘，
{A0A1A2A3}={11111−i−1i1−11−11i−1−i}{a0a1a2a3}\begin{Bmatrix} A_0 \\ A_1 \\ A_2 \\ A_3 \end{Bmatrix} = \begin{Bmatrix} 1 & 1 & 1 & 1 \\ 1 & -i & -1 & i \\ 1 & -1 & 1 & -1 \\ 1 & i & -1 & -i \end{Bmatrix} \begin{Bmatrix} a_0 \\ a_1 \\ a_2 \\ a_3 \end{Bmatrix} ⎩⎪⎪⎨⎪⎪⎧A0A1A2A3⎭⎪⎪⎬⎪⎪⎫=⎩⎪⎪⎨⎪⎪⎧11111−i−1i1−11−11i−1−i⎭⎪⎪⎬⎪⎪⎫⎩⎪⎪⎨⎪⎪⎧a0a1a2a3⎭⎪⎪⎬⎪⎪⎫

为了快速计算AAA，我们可以将计算稍作改变：
A0=(a0+a2)+(a1+a3)A1=(a0−a2)−i(a1−ia3)A2=(a0+a2)−(a1+a3)A3=(a0−a2)+i(a1−a3)A_0 = (a_0+a_2) + (a_1+a_3) \newline A_1 = (a_0-a_2) - i(a_1-ia_3) \newline A_2 = (a_0+a_2) - (a_1+a_3) \newline A_3 = (a_0-a_2) + i(a_1-a_3) A0=(a0+a2)+(a1+a3)A1=(a0−a2)−i(a1−ia3)A2=(a0+a2)−(a1+a3)A3=(a0−a2)+i(a1−a3)
这样很多运算结果可以公用。

如果使用图线的方式来表示乘与加，有如下表示方法，

所以4个sample的DFT可以画成如下，

这便是DFT中蝴蝶算法的基本原理。
下面是8个sample的运算图

自己实现DFT

根据上述算法，在matlab上做如下DFT实现

定义DFT函数

function [f, Xk] = mydft(xn, fs, N)
%UNTITLED2 此处显示有关此函数的摘要
%   此处显示详细说明
%   xn:  the discrete samples
%   fs:  frequency
%   N:   the amount of samplesn = 0:1:N-1;k = n;WN= exp(-1i*2*pi/N);nk = n' * k;Xk = xn(1:N) * WN.^nk;f = 0:fs/N:fs/N*(N-1);
end

Demo代码如下，

可以看到和Matlab内置函数fft结果是一样的

复数域的三维观测

频率在时间轴上形成了“波”形的二维形状，但是它仅仅展示出了振幅（能量）信息。
傅里叶变化增加了复数域轴，使它展示出了相位信息，每一个点也成为了三维空间上的点。

下面我们在三维视角下观测，操作步骤如下，

绘制波形

N=200;fs=20;
i=0:1:N-1;
xn = sin(2*pi*i*fs/N);

绘制振幅谱

Xk = fft(xn,N);
mag = abs(Xk);plot(1:1:N/2, mag(1:N/2));

绘制real部分

scatter(1:1:N, real(Xk));

绘制img部分

scatter(1:1:N, imag(Xk));

投影图

scatter(real(Xk), imag(Xk));

绘制3D点

scatter3(1:1:N, real(Xk), imag(Xk));

有一些有意思的结论，

real部分是对称的，img部分是共轭兑成的
如果是完整的sin和cos波，在频率点处，可以看到跳点，原来的“直线”像断在这里一样
如果是完整的sin和cos波，投影图是趋向于X型的

选取一段钢琴音，观测它的3D图

real部分图

imagine部分图

投影图基本呈现X形状

附录

8sample DFT 样例

旋转因子

W81=e−i2π/8=e−iπ/4=22−i22W_8^1 = e^{-i2\pi/8} = e^{-i\pi/4} = \frac{\sqrt{2}}2-i\frac{\sqrt{2}}2W81=e−i2π/8=e−iπ/4=22−i22
W82=−iW_8^2 = -iW82=−i
W83=−22−i22W_8^3 = -\frac{\sqrt{2}}2-i\frac{\sqrt{2}}2W83=−22−i22
W84=−1W_8^4 = -1W84=−1
W85=−22+i22W_8^5 = -\frac{\sqrt{2}}2+i\frac{\sqrt{2}}2W85=−22+i22
W86=iW_8^6 = iW86=i
W87=22+i22W_8^7 = \frac{\sqrt{2}}2+i\frac{\sqrt{2}}2W87=22+i22
W88=1W_8^8 = 1W88=1

旋转因子矩阵

{W80W80W80W80W80W80W80W80W80W81W82W83W84W85W86W87W80W82W84W86W80W82W84W86W80W83W86W81W84W87W82W85W80W84W80W84W80W84W80W84W80W85W82W87W84W81W86W83W80W86W84W82W80W86W84W82W80W87W86W85W84W83W82W81}\begin{Bmatrix} W_8^0\space\space W_8^0\space\space W_8^0\space\space W_8^0\space\space W_8^0\space\space W_8^0\space\space W_8^0\space\space W_8^0 \\ \\ W_8^0\space\space W_8^1\space\space W_8^2\space\space W_8^3\space\space W_8^4\space\space W_8^5\space\space W_8^6\space\space W_8^7 \\ \\ W_8^0\space\space W_8^2\space\space W_8^4\space\space W_8^6\space\space W_8^0\space\space W_8^2\space\space W_8^4\space\space W_8^6 \\ \\ W_8^0\space\space W_8^3\space\space W_8^6\space\space W_8^1\space\space W_8^4\space\space W_8^7\space\space W_8^2\space\space W_8^5 \\ \\ W_8^0\space\space W_8^4\space\space W_8^0\space\space W_8^4\space\space W_8^0\space\space W_8^4\space\space W_8^0\space\space W_8^4 \\ \\ W_8^0\space\space W_8^5\space\space W_8^2\space\space W_8^7\space\space W_8^4\space\space W_8^1\space\space W_8^6\space\space W_8^3 \\ \\ W_8^0\space\space W_8^6\space\space W_8^4\space\space W_8^2\space\space W_8^0\space\space W_8^6\space\space W_8^4\space\space W_8^2 \\ \\ W_8^0\space\space W_8^7\space\space W_8^6\space\space W_8^5\space\space W_8^4\space\space W_8^3\space\space W_8^2\space\space W_8^1 \\ \\ \end{Bmatrix} ⎩⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎨⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎧W80 W80 W80 W80 W80 W80 W80 W80W80 W81 W82 W83 W84 W85 W86 W87W80 W82 W84 W86 W80 W82 W84 W86W80 W83 W86 W81 W84 W87 W82 W85W80 W84 W80 W84 W80 W84 W80 W84W80 W85 W82 W87 W84 W81 W86 W83W80 W86 W84 W82 W80 W86 W84 W82W80 W87 W86 W85 W84 W83 W82 W81⎭⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎬⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎫

Discrete Fourier Transform离散傅里叶变换算法相关推荐

理解快速离散傅里叶变换算法(FFT)
本文是视频The Fast Fourier Transform (FFT): Most Ingenious Algorithm Ever?的整理. 离散傅里叶变换的用法 FFT是一个非常快速的离散傅里 ...
Fast Fourier transform快速傅里叶变换
0 傅里叶分析和滤波资料来源:https://ww2.mathworks.cn/help/matlab/fourier-analysis-and-filtering.html?s_tid=CRUX_ ...
基于CUDA的离散傅里叶变换（Discrete Fourier Transform，DFT）
最近在做地震勘探的全波形反演,用分频反演的方法,需要对地震波场按照特定的频段进行傅里叶变换,这要用到DFT.在CPU上,DFT的计算非常耗时,当处理三维数据时耗时更加严重,所以,本人用CUDA+SU( ...
Discrete Fourier Transform（DCT）的理解
参考链接:简要理解DFT_我叫夏满满的博客-CSDN博客_dft原理 FFT(快速傅里叶变换)其本质就是DFT,只不过可以快速的计算出DFT结果.要理解FFT,必须先理解DFT(DiscreteFou ...
SAR ADC系列2：DFT离散傅里叶变换
目录 ADC动态性能仿真/测试平台 DFT:离散傅里叶分析 DFT的 Matlab 实现: 频谱泄露: 如何规避频谱泄露: ADC性能分析:DFT Cadence环境下的DFT分析实例: Matlab ...
傅里叶变换一维快速傅里叶变换（快速的一维离散傅里叶变换、分治法）
一.介绍 1.一维离散傅里叶变换DFT. DFT:(Discrete Fourier Transform)离散傅里叶变换是傅里叶变换在时域和频域上都呈离散的形式,将信号的时域采样变换为其DTFT的频域 ...
傅里叶变换一维离散傅里叶变换
1.介绍. DFT:(Discrete Fourier Transform)离散傅里叶变换是傅里叶变换在时域和频域上都呈离散的形式,将信号的时域采样变换为其DTFT的频域采样.在形式上,变换两端(时域 ...
傅里叶变换二维离散傅里叶变换
1.介绍. DFT:(Discrete Fourier Transform)离散傅里叶变换是傅里叶变换在时域和频域上都呈离散的形式,将信号的时域采样变换为其DTFT的频域采样.在形式上,变换两端(时域 ...
Python+OpenCV：傅里叶变换(Fourier Transform)
Python+OpenCV:傅里叶变换(Fourier Transform) ############################################################# ...

Discrete Fourier Transform离散傅里叶变换算法

DFT

iDFT

具体示例

2个sample的DFT

4个sample的DFT

自己实现DFT

定义DFT函数

复数域的三维观测

绘制波形

绘制振幅谱

绘制real部分

绘制img部分

投影图

绘制3D点

附录

8sample DFT 样例

旋转因子矩阵

Discrete Fourier Transform离散傅里叶变换算法相关推荐

最新文章

热门文章