01 - 关于傅里叶变换及 FFT

CliCli_0

于 2024-03-22 11:33:19 发布

阅读量977

点赞数 19

分类专栏：笔记本文章标签：傅里叶分析

本文链接：https://blog.csdn.net/CliCli_0/article/details/136935787

版权

笔记本专栏收录该内容

5 篇文章 0 订阅

订阅专栏

1. 傅里叶变换说明

1.1. 中英文对照表

英文简称	英文全称	说明
FT	Fourier Transformation	傅里叶变换，针对连续非周期信号的傅里叶分析
FS	Fourier Serie	傅里叶级数，任意连续周期信号在各正弦信号分量上的权重
DTFT	Discrete-time Fourier Transform	离散时间傅里叶变换，离散非周期信号的 FT
DFT	Discrete Fourier Transform	离散傅里叶变换，离散周期信号的 FT
FFT	Fast Fourier Transformation	快速傅里叶变换

1.2 傅里叶相关变换公式

FT （连续非周期信号）：

$F(\omega)= \int_{-\infty}^{\infty}f(t)e^{-j{\omega}t}dt$

DTFT （连续非周期信号进行采样，然后做 FT ）：

$X(\omega) = \sum_{n=-\infty}^{+\infty}x(n)e^{-j\omega n}$

DFT （在频域区间 [0, 2 $\pi$ ] 内以 $\frac{2\pi}N$ 为间隔对 DTFT 进行采样）：

$\sum_{n=0}^{N-1}x(n)e^{-j\frac{2\pi}Nkn}$

IDFT （DFT 反变换，注意系数 $\frac1N$ ）
$\frac1N\sum_{k=0}^{N-1}X(k)e^{j\frac{2\pi}{N}nk}$
FFT
- 相关性质
  
  令 $W_N^{kn} = e^{-j\frac{2\pi}Nkn}$ ，则：
  
  1). $W_N^0 = W_N^N = 1$
  
  2). $W_N^k = W_N^{k+N}$
  
  3). $W_{dN}^{dk} = W_{N}^{k}$
  
  4). $W_2^1 = -1$
  
  5). $W_N^{k+\frac{N}2} = -W_N^k$
  
  6). $W_N^{k+n} = W_N^k * W_N^n$
- 公式推导
  
  设有如下多项式，其中 $n$ 为 2 的整数次幂：
  $f(x) = a_0 + a_1x^1 + a_2x^2 + a_3x^3 + ... + a_{n-1}x^{n-1}$
  分离偶数次幂项与奇数次幂项，并令：
  $f_1(x) = a_0 + a_2x^2 + a_4x^4 + ... + a_{n-2}x^{n-2}$
  
  $f_2(x) = a_1 + a_3x^2 + a_5x^4 + ... + a_{n-1}x^{n-2}$
  
  则有：
  $f(x) = f'(x^2)= f_1(x^2) + xf_2(x^2)$
  令 $x = w_n^k$ ，且 $\frac{n}2$ 时，由性质 3 可得
  $f(w_n^k) = f'((w_n^k)^2) = f_1((w_n^k)^2) + w_n^kf_2((w_n^k)^2) = f_1((w_{n/2}^k)) + w_n^kf_2((w_{n/2}^k))$
  当 $k^{'}$ 取后半部分，即 $\frac{n}2$ 时，由性质 2 、 5 得
  $f(w_n^{k'}) = f_1((w_n^{k+{n/2}})^2) + w_n^{k+{n/2}}f_2((w_n^{k+{n/2}})^2) = f_1((w_{n/2}^k)) - w_n^kf_2((w_{n/2}^k))$
  结合 DFT 公式，可知
  
  1）n 点的 DFT 可以转换成 n/2 的 DFT
  
  2）计算前一半 DFT 值时，即可同时得到后一半 DFT 值
- 矩阵示例
  
  对信号 $x (t)$ 在一个周期内进行 $N$ 次采样，则 $N$ 点的 DFT 满足：
  
  $\begin{Bmatrix} w_n^{00}&w_n^{01}&\cdots&w_n^{0(n-1)}\\ w_n^{10}&w_n^{11}&\cdots&w_n^{1(n-1)}\\ \vdots&\vdots&\ddots&\vdots\\ w_n^{(n-1)0}&w_n^{(n-1)1}&\cdots&w_n^{(n-1)(n-1)} \end{Bmatrix} * \begin{Bmatrix} x_0\\ x_1\\ \vdots\\ x_{n-1} \end{Bmatrix} = \begin{Bmatrix} X_0\\ X_1\\ \vdots\\ X_{n-1} \end{Bmatrix}$
  
  $n$ 表示采样点下标，即 $w$ 矩阵列数； $k$ 表示周期内的频率分量下标，即 $w$ 矩阵行数
  
  将 $w$ 矩阵偶、奇列重排，DFT 等式可表示为
  
  $\begin{Bmatrix} x_0\\ x_2\\ \vdots\\ x_{n-2}\\ x_1\\ x_3\\ \vdots\\ x_{n-1} \end{Bmatrix} = \begin{Bmatrix} X_0\\ X_2\\ \vdots\\ X_{n-2}\\ X_1\\ X_3\\ \vdots\\ X_{n-1} \end{Bmatrix}$
  
  此时，
  $\begin{Bmatrix} w_n^{00}&w_n^{02}&\cdots&w_n^{0(n-2)}&w_n^{01}&w_n^{03}&\cdots&w_n^{0(n-1)}\\ w_n^{10}&w_n^{12}&\cdots&w_n^{1(n-2)}&w_n^{11}&w_n^{13}&\cdots&w_n^{1(n-1)}\\ \vdots&\vdots&\ddots&\vdots&\vdots&\vdots&\ddots&\vdots\\ w_n^{(\frac{n}2-1)0}&w_n^{(\frac{n}2-1)2}&\cdots&w_n^{(\frac{n}2-1)(n-2)}&w_n^{(\frac{n}2-1)1}&w_n^{(\frac{n}2-1)3}&\cdots&w_n^{(\frac{n}2-1)(n-1)}\\ w_n^{\frac{n}20}&w_n^{\frac{n}22}&\cdots&w_n^{\frac{n}2(n-2)}&w_n^{\frac{n}21}&w_n^{\frac{n}23}&\cdots&w_n^{\frac{n}2(n-1)}\\ w_n^{(\frac{n}2+1)0}&w_n^{(\frac{n}2+1)2}&\cdots&w_n^{(\frac{n}2+1)(n-2)}&w_n^{(\frac{n}2+1)1}&w_n^{(\frac{n}2+1)3}&\cdots&w_n^{(\frac{n}2+1)(n-1)}\\ \vdots&\vdots&\ddots&\vdots&\vdots&\vdots&\ddots&\vdots\\ w_n^{(n-1)0}&w_n^{(n-1)2}&\cdots&w_n^{(n-1)(n-2)}&w_n^{(n-1)1}&w_n^{(n-1)3}&\cdots&w_n^{(n-1)(n-1)} \end{Bmatrix}$
  
  将 $W^{'}$ 矩阵划分为四部分
  $\begin{Bmatrix} A&B\\ C&D \end{Bmatrix}$
  
  其中， $A$ 矩阵中各项为
  $\begin{Bmatrix} w_n^{00}&w_n^{02}&\cdots&w_n^{0(n-2)}\\ w_n^{10}&w_n^{12}&\cdots&w_n^{1(n-2)}\\ \vdots&\vdots&\ddots&\vdots\\ w_n^{(\frac{n}2-1)0}&w_n^{(\frac{n}2-1)2}&\cdots&w_n^{(\frac{n}2-1)(n-2)} \end{Bmatrix}$
  
  则 1）依据性质 2 可得
  $C = A$
  
  2）依据性质 6 可得（注意此处 $k < n /2$ ）
  
  $\begin{Bmatrix} w_n^{0}&\cdots&\cdots&\cdots\\ \cdots&w_n^{1}&\cdots&\cdots\\ \cdots&\cdots&w_n^{k}&\cdots\\ \cdots&\cdots&\cdots&w_n^{(\frac{n}2-1)}\\ \end{Bmatrix} *A = diag(w_n^k)*A$
  
  3）依据性质 2 5 可得
  $D = -B = -diag(w_n^k)*A$
  
  4）依据性质 3 可将 $A$ 矩阵简化为
  $\begin{Bmatrix} w_{\frac{n}2}^{00}&w_{\frac{n}2}^{01}&\cdots&w_{\frac{n}2}^{0(\frac{n}2-1)}\\ w_{\frac{n}2}^{10}&w_{\frac{n}2}^{11}&\cdots&w_{\frac{n}2}^{1(\frac{n}2-1)}\\ \vdots&\vdots&\ddots&\vdots\\ w_{\frac{n}2}^{(\frac{n}2-1)0}&w_{\frac{n}2}^{(\frac{n}2-1)1}&\cdots&w_{\frac{n}2}^{(\frac{n}2-1)(\frac{n}2-1)} \end{Bmatrix}$
  此时 $A$ 即 $n /2$ 点 DFT 的 $W_{\frac{n}2}$ ，对于 $n$ 点的 DFT 可以拆分为 $n /2$ 点的 DFT ，且系数因子（即旋转因子）满足：
  $W_n = \begin{Bmatrix} W_{\frac{n}2}&diag(w_n^k)*W_{\frac{n}2}\\ W_{\frac{n}2}&-diag(w_n^k)*W_{\frac{n}2} \end{Bmatrix}$
  
  特别地，对于 2 点的 DFT 有
  $\begin{Bmatrix} X_0\\ X_1 \end{Bmatrix} = W_2 * \begin{Bmatrix} x_0\\ x_1 \end{Bmatrix} = \begin{Bmatrix} 1&1\\ 1&-1 \end{Bmatrix} * \begin{Bmatrix} x_0\\ x_1 \end{Bmatrix}$
  
  8 点 DFT 蝶形运算图示
- 代码实现
```
function y=myditfft (x)
%本程序对输入序列 x实现DIT-FFT基2算法，点数取大于等于x长度的2的幂次
%y=myditfft (x)
x=randn(1,2^4);
m=nextpow2 (length(x)); %进行补零操作，满足2的整数幂----nextpow2（7）=3，
N=2 ^m;         %求x的长度对应的2的最低幂次m，对应多少点的FFT

if length(x) <N
    x=[x, zeros(1, N-length(x))] ;%若x的长度不是2的幂，补零到2的整数幂.
end
nxd=bin2dec (fliplr (dec2bin((1:N)-1,m)))+1;    %求1:2^m数列的倒序
%bin2dec将二进制转换为10进制   flip1r倒叙排列   dec2bin把一个十进制数转换成一个字符串形式的二进制数
y=x(nxd) ;          %将x倒序排列作为y的初始值，这样我们输出的数值和蝶形运算的数值相对应。

for mm=1 :m         %将DFT作m次基2分解,从左到右，对每次分解作DFT运算——把DFT分成log2N段来实现
    Nmr= 2 ^mm;     %这里是旋转因子中的N    
    u=1;            %初始化WN0=1
    WN=exp (-1i*2*pi/Nmr) ;     %本次分解的基本DFT因子WN=exp (-i*2*pi/ Nmr)
    for j=1 :Nmr/2              %本次跨越间隔内的各次蝶形运算
        for k=j :Nmr :N         %本次蝶形运算的跨越间隔为Nmr=2 ^mm
            kp=k+Nmr/2;         %确定蝶形运算的对应单元下标.
            t=y (kp) *u;        %蝶形运算的乘积项
            y (kp)=y(k) -t;     %蝶形运算的减法项
            y(k)=y(k) +t;       %蝶形运算的加法项
        end
        u=u*WN;                 %修改旋转因子，多乘一个基本DFT因子WN
    end
end
```
  【参考1】： FFT 代码实现
  
  【参考2】：蝶形运算输入数据重排原理
  
  【参考3】：图示蝶形运算