线性变换+DFT+滤波
最近又点开了一个关于傅立叶变换的文章，里面通过动画的形式展示了如何将一个时域的输入信号展开成多个正余弦信号的叠加。看起来好像醍醐灌顶，懂了，然后又忘了。
其实，仔细想了想，傅立叶变换变的是坐标系。我们拿离散数据来讲，时域的输入信号其实是一个处于标准坐标系内的（列）向量 x ∈ R N × 1 \mathbf{x} \in \mathbb{R}^{N\times 1} x∈RN×1， 而傅立叶变换就是通过线性变换的方式改变坐标系，形成一个新的向量： x ^ = F x \hat{\mathbf{x}} = \mathbf{F}\mathbf{x} x^=Fx。
让我们先回顾一下矩阵线性变换吧，具体可以看之前的文章：矩阵于线性代数。
线性变换
假设我们有一个向量
v
=
[
1
;
2
]
\mathbf{v} = [1; 2]
v=[1;2]，通常，1和2表示的是这个向量在标准坐标系两个轴上的投影，而这个标准的坐标系也是可以用矩阵，
M
\mathbf{M}
M，来表示：
M
=
[
1
0
0
1
]
\mathbf{M} = \begin{bmatrix}1 & 0\\ 0 & 1\end{bmatrix}
M=[1001]。如果这个向量是基于另一个坐标系的，比方说：
[
0
−
1
1
0
]
\begin{bmatrix}0 & -1\\ 1 & 0\end{bmatrix}
[01−10] （坐标系逆时针旋转90度），我们可以通过以下计算，求出这个向量在标准坐标系下的投影：
v
′
=
[
0
1
−
1
0
]
v
=
[
−
2
1
]
\mathbf{v}^{'} = \begin{bmatrix}0 & 1\\ -1 & 0\end{bmatrix}\mathbf{v} = \begin{bmatrix}-2 \\ 1\end{bmatrix}
v′=[0−110]v=[−21]
线性变换的意义是将一个在标准坐标系中不太容易分析的数据放在另一个坐标系下，而这另一个坐标系往往可以让数据分析变得特别简单。
DFT+滤波
回到离散傅立叶变换，标准坐标系会被转换成新的坐标系，由矩阵
F
\mathbf{F}
F 表示。新的坐标系下的向量
x
\mathbf{x}
x 又被投影到标准坐标系中，得到了
x
^
\hat{\mathbf{x}}
x^。至于矩阵
F
\mathbf{F}
F，所有相关的书中都是有它的表达式的：
[
1
1
1
⋯
1
1
ω
n
ω
n
2
⋯
ω
n
n
−
1
1
ω
n
2
ω
n
4
⋯
ω
n
2
(
n
−
1
)
⋮
⋮
⋮
⋱
⋮
1
ω
n
n
−
1
ω
n
2
(
n
−
1
)
⋯
ω
n
(
n
−
1
)
2
]
\begin{bmatrix} 1 & 1 & 1 & \cdots & 1\\ 1 & \omega_n & \omega_n^2 & \cdots & \omega_n^{n-1}\\ 1 & \omega_n^2 & \omega_n^4 & \cdots & \omega_n^{2(n-1)}\\ \vdots & \vdots & \vdots & \ddots & \vdots\\ 1 & \omega_n^{n-1} & \omega_n^{2(n-1)} & \cdots & \omega_n^{(n-1)^2} \end{bmatrix}
⎣⎢⎢⎢⎢⎢⎢⎡111⋮11ωnωn2⋮ωnn−11ωn2ωn4⋮ωn2(n−1)⋯⋯⋯⋱⋯1ωnn−1ωn2(n−1)⋮ωn(n−1)2⎦⎥⎥⎥⎥⎥⎥⎤
其中，
ω
n
=
e
−
2
π
i
/
n
\omega_n=e^{-2\pi i/n}
ωn=e−2πi/n。至于如何从正余弦函数的表示通过欧拉变换，得到这个
ω
\omega
ω 就不再赘述了。
在MATLAB计算的时候，比方说，我们随机生成一个长度为 50 50 50 的时域信号 s \mathbf{s} s，我们可以通过上面的表达式以及已知的信号长度来具体表示矩阵 F \mathbf{F} F，然后就可以求出 s ^ \hat{\mathbf{s}} s^ 了：
s = rand(50, 1);
n = length(s);
w = exp(-1i*2*pi/n);
[I,J] = meshgrid(1:n,1:n);
F = w.^((I-1).*(J-1));
s_h = F*s;
当然，这种计算的方式在实际应用时往往耗时，所以有了快速傅立叶变换！下面给出一个使用FFT和基于DFT的例子：
Fs = 96e3; % Sampling frequency
T = 1/Fs; % Sampling period
L = 1024; % Length of signal
t = (0:L-1)*T; % Time vector
S = 0.7*sin(2*pi*5e3*t) + sin(2*pi*12e3*t) + sin(2*pi*18e3*t) + sin(2*pi*26e3*t);
subplot(3,1,1), plot(Fs*t(1:200), S(1:200))
f = Fs*(0:(L/2))/L;
tic
Y = fft(S);
timeElapsed = toc
P2 = abs(Y/L);
P1 = P2(1:L/2+1);
P1(2:end-1) = 2*P1(2:end-1);
subplot(3,1,2), plot(f,P1)
tic
n = length(S);
w = exp(-1i*2*pi/n);
[I,J] = meshgrid(1:n,1:n);
DFT = w.^((I-1).*(J-1));
Y2 = DFT*S';
timeElapsed = toc
P2 = abs(Y2/L);
P1 = P2(1:L/2+1);
P1(2:end-1) = 2*P1(2:end-1);
subplot(3,1,3), plot(f,P1)
通过MATLAB的tic’toc来记录运算时间，我的笔记本给出的时间是0.0032和0.1162秒，对应于FFT和DFT。代码中的注意点：
- 采样频率 F s F_s Fs 和样本长度 L L L 决定了每条谱线的分辨率。在上面的例子中，傅立叶变换后得到的数据长度也是 L L L，其中，每个点之间间隔： F s / L = 93.75 F_s/L = 93.75 Fs/L=93.75 Hz。
- 能显示的频谱条数是所选择的傅立叶变换的输入长度的一半（采样定律）。
基于这些知识，我们可以很容易在实践中运用傅立叶变换。比方说，我们可以通过示波器的FFT功能来观测某个信号，通过频谱图，我们可以通过修改电路，过滤某些噪音信号。如果信号受到白噪音的干扰，比方如下：
Fs = 20e3; % Sampling frequency
T = 1/Fs; % Sampling period
L = 1024; % Length of signal
t = (0:L-1)*T; % Time vector
S = sin(2*pi*1e3*t) + sin(2*pi*5e3*t);
S_n = S + 2.5*randn(size(t)); % Add some noise
subplot(3,1,1)
plot(Fs*t(1:200), S(1:200),'k','LineWidth',1.2), hold on
plot(Fs*t(1:200), S_n(1:200),'r','LineWidth',1.5)
legend('Clean','Noisy')
f = Fs*(0:(L/2))/L;
Y = fft(S_n);
PSD = Y.*conj(Y)/L; % Power spectrum (power per freq)
%Use the PSD to filter out noise
indices = PSD>100; % Find all freqs with large power
PSDclean = PSD.*indices; % Zero out all others
subplot(3,1,2)
plot(f,PSD(1:L/2+1),'r','LineWidth',1.5), hold on
plot(f,PSDclean(1:L/2+1),'b','LineWidth',1.2)
legend('Noisy','Filtered')
Y = indices.*Y;
S_f = ifft(Y); % Inverse FFT for filtered time signal
subplot(3,1,3)
plot(Fs*t(1:200), S(1:200),'k','LineWidth',1.2), hold on
plot(Fs*t(1:200), S_f(1:200),'b','LineWidth',1.2)
legend('Clean','Filtered')
运行以上代码可以得到下图：
我们可以看到，加入高斯噪音后，原本信号杂乱无章。通过傅立叶变换，我们看到除了两个频率点有强功率外，其余都是很弱的信号。我们可以通过设定一个阀值来过滤那些弱信号。逆变换之后得到的信号基本与原本的信号一致。