MATLAB图解傅里叶变换（初学者也可以理解）

1、概述

相信很多人对于傅里叶变换可能觉得比较复杂和有点难懂，其实不难，它只是一种积分变换。
傅里叶变换，表示能将满足一定条件的某个函数表示成三角函数（正弦和/或余弦函数）或者它们的积分的线性组合。也就是说"任意"的函数通过一定的分解，都能够表示为正弦函数的线性组合的形式。而正弦函数在物理上是被充分研究而相对简单的函数类。因为特别好使，所以傅里叶变换在物理学、数论、组合数学、信号处理、概率、统计、密码学、声学、光学等领域都有着广泛的应用。
为什么使用正弦曲线来分解原函数呢？
因为正弦曲线的保真度。一个正弦曲线信号输入后，输出的仍是正弦曲线，只有幅度和相位可能发生变化，但是频率和波的形状仍是一样的。且只有正弦曲线才拥有这样的性质，所以我们会将原函数分解成正弦函数来处理。

2、连续傅里叶变换Continuous Fourier Transform,CFT

平时看到的傅里叶变换如果没有前缀，通常指的就是连续傅里叶变换。我们先来看一张图片：

从图片中我们可以了解到傅里叶的变换是时间域和频率域的互相转换，由于是连续的，所以使用的是积分。
平时我们所听到的音乐，就是随着时间的变化，听到不同的声音，乐器声和人声的各种混合，如果我想单独将人声给提取出来，这个时候我们可以使用傅里叶变换来处理，将这些不同的声音分解成各自的频率所表示的形式。这样我们就将本来很难处理的时域信号转换成了容易处理的频域信号了（信号的频谱），再利用一些工具对这些频域信号进行处理、加工即可。最后还可以利用傅里叶反变换将这些频域信号转换成时域信号。

3、离散傅里叶变换Discrete Fourier Transform,DFT

实际上我们的计算机是只能处理离散的数值信号，跟数学上的计算是有区别的，这些听到的连续的声音就需要做离散处理，因为计算机本质上只能处理0和1，所以这里不能使用积分，我们知道积分其实是极限的连续求和对吧，那我们这里就只能是离散求和了，而且不能处理无穷，计算机只能处理有限性，这个是前提条件。我们来看一张图片：

可以看到我们的离散傅里叶使用的是求和符号而不是上面的积分符号，在计算机里面我们就是使用这个DFT来处理信号。
直接使用DFT的定义计算的计算复杂度为 O(N²)，而快速傅里叶变换(Fast Fourier Transform, FFT)可以将复杂度改进为 O(nlogn)。计算复杂度的降低以及数字电路计算能力的发展，使得DFT成为在信号处理领域十分实用且重要的方法。

从上面的分析来看，我们知道时域分析与频域分析是对信号的两个观察面。时域分析是以时间轴为坐标表示动态信号的关系；频域分析是把信号变为以频率轴为坐标表示出来。
一般来说，时域的表示较为形象与直观，频域分析则更为简练，剖析问题更为深刻和方便。信号分析的趋势是从时域向频域发展。然而，它们是互相联系，缺一不可，相辅相成的，它们的转换就是我们讲的傅里叶变换。

4、MATLAB应用

我们来看下具体的应用，在MATLAB中使用快速傅里叶变换FFT，将会是怎么样的。

先来看下fft函数的定义：

Y = fft(X) 用快速傅里叶变换 (fft) 算法计算 X 的离散傅里叶变换 (DFT)。
如果 X 是向量，则 fft(X) 返回该向量的傅里叶变换。
如果 X 是矩阵，则 fft(X) 将 X 的各列视为向量，并返回每列的傅里叶变换。
如果 X 是一个多维数组，则 fft(X) 将沿大小不等于 1 的第一个数组维度的值视为向量，并返回每个向量的傅里叶变换。

Y = fft(X,n) 返回 n 点 DFT。如果未指定任何值，则 Y 的大小与 X 相同。
如果 X 是向量且 X 的长度小于 n，则为 X 补上尾零以达到长度 n。
如果 X 是向量且 X 的长度大于 n，则对 X 进行截断以达到长度 n。
如果 X 是矩阵，则每列的处理与在向量情况下相同。
如果 X 为多维数组，则大小不等于 1 的第一个数组维度的处理与在向量情况下相同

Y = fft(X,n,dim) 返回沿维度 dim 的傅里叶变换。X 是矩阵，比如：
fft(X,n,1) 沿 X 的各列进行运算，返回每列的傅里叶变换；
fft(X,n,2) 沿 X 的各行进行运算，返回每行的傅里叶变换。

4.1、噪声信号Noise signal

我们来看具体的示例，使用傅里叶变换求噪声中隐藏的信号的频率分量（频谱分析）。

Fs = 1000; %采样频率
T = 1/Fs; %采样周期
L = 1500; %信号长度(1.5秒)
t = (0:L-1)*T; %时间

%创建一个信号，其中包含幅值为0.7，频率为50Hz和幅值为1，频率为120Hz的正弦波
S = 0.7*sin(2*pi*50*t) + sin(2*pi*120*t);
%添加一个标准差为2的噪声
X = S + 2*randn(size(t));
%绘制噪声时域信号。
plot(1000*t(1:50),X(1:50))
title('加了噪声的信号')
xlabel('t (毫秒)')
ylabel('X(t)')

通过查看这个时域信号X(t)，我们很难确定频率分量，所以我们通过计算傅里叶变换来看频谱分析，这里使用的是快速傅里叶变换： Y = fft(X);

%计算双侧频谱P2，然后基于P2和偶数信号长度L计算单侧频谱P1，单侧频谱表示信号频率的幅度分布只从0到正频率一侧，正负频率是对称的
P2 = abs(Y/L);
P1 = P2(1:L/2+1);
P1(2:end-1) = 2*P1(2:end-1);

其中这里的abs是返回绝对值，但是这里是复数的形式，所以返回的是它们的模长(幅值)，表示的是复数的大小，对于复数还有一个属性比较重要就是相位，表示的是复数的方向。
我们画图看下频谱对应的幅值：

f = Fs*(0:(L/2))/L;
plot(f,P1) 
title('单侧幅值频谱')
xlabel('f (Hz)')
ylabel('|P1(f)|')

可以看到在50Hz和120Hz的地方，幅值基本上还是接近0.7和1的，因为加了噪声的缘故，如果没有噪声那肯定是对应着0.7和1，这个大家可以试下。这里我们就可以看出，在时域坐标中，查看信号的频率分量比较困难，当我们使用傅里叶变换之后，在频域中查看信号就变得简单了。

4.2、高斯脉冲Gaussian Pulses

所谓的脉冲就是发送信号，接通一下，断开一下，这个就是一个脉冲，现在的步进电机和伺服电机的驱动都是通过脉冲来驱动。
高斯脉冲是具有高斯形状的时间强度曲线的脉冲，其中高斯函数如下：

$f(x) = a*e^{\frac{-(x - b)^2}{2c^2}}$

其中a、b、c是常数，a决定了曲线的高度(峰值)，b决定了曲线的中心位置(横坐标)，c决定了曲线的宽度(标准差)。
因为高斯脉冲的时间强度曲线呈现出高斯分布的特点，所以它在时间和频率上都能达到较好的性能，从而在许多领域得到广泛的应用。
一起来看下高斯脉冲在时域转到频域的情况：

Fs = 100; %采样频率
t = -0.5:1/Fs:0.5; %时间
L = length(t); %信号长度
%高斯脉冲对应上面公式中的a、b、c了解
X = 1/(4*sqrt(2*pi*0.01))*(exp(-t.^2/(2*0.01)));

%画图
plot(t,X)
title('时域中的高斯脉冲')
xlabel('Time (t)')
ylabel('X(t)')

从图形中，我们可以看到是一个高斯分布，其中心点是0，峰值约为0.9974，另外公式看出宽度是0.1，这个函数在傅里叶变换中属于典型的特征函数。
在转换频域之前，首先从原始信号长度确定是下一个2次幂的新输入长度。当信号长度不等于2的幂时，使用 nextpow2 函数可提升 fft 的性能，使用的是尾随零填充信号。

n = 2^nextpow2(L);
%将高斯脉冲转换为频域。
Y = fft(X,n);%length(X):101，n:128，所以length(Y)是128
%画图
f = Fs*(0:(n/2))/n;
P = abs(Y/n);

plot(f,P(1:n/2+1)) 
title('频域中的高斯脉冲')
xlabel('Frequency (f)')
ylabel('|P(f)|')

代码中的nextpow2函数，返回的是下一个2次幂大于等于这些值(绝对值)的指数，比如：

a = [1 -2 3 -4 5 9 519];
p = nextpow2(a)
p =
     0     1     2     2     3     4    10

2.^p 或者 pow2(p)
ans =
     1      2      4      4      8     16    1024

4.3、余弦波

我们来看下最常见的一个波形，比较它们不同频率的情况：

Fs = 1000; %采样频率
T = 1/Fs; %采样周期
L = 1000; % 信号长度
t = (0:L-1)*T; % 时间

%三个余弦波，区别就是频率不一样
x1 = cos(2*pi*50*t); %50Hz
x2 = cos(2*pi*150*t); %150Hz
x3 = cos(2*pi*300*t); %300Hz
% 组合成一个矩阵
X = [x1; x2; x3];

%画出每个余弦波前100个条目
for i = 1:3
    subplot(3,1,i)
    plot(t(1:100),X(i,1:100))
    title(['余弦波 ',num2str(i),' 的时域图形'])
end

这里我们发现三个余弦波，只能说往下的余弦波的频率比上面的要大，还是不能很清晰的表现出其本质，我们将其转换到频域空间看下是什么效果：

%同样的性能考虑，使用nextpow2函数定义新长度
n = 2^nextpow2(L);
%接下来就计算傅里叶变换
dim = 2;
Y = fft(X,n,dim);

%计算每个信号的双侧频谱和单侧频谱
P2 = abs(Y/L);
P1 = P2(:,1:n/2+1);
P1(:,2:end-1) = 2*P1(:,2:end-1);

%在频域图形里面，画出每一个余弦波的单侧幅值频谱
for i=1:3
    subplot(3,1,i)
    plot(0:(Fs/n):(Fs/2-Fs/n),P1(i,1:n/2))
    title(['余弦波 ',num2str(i),'在频域中的图形'])
end