Matlab可以说是一个非常有用且功能齐全的工具，在通信、自控、金融等方面有广泛的应用。
本文讨论使用Matlab对信号进行频域分析的方法。
说到频域，不可避免的会提到傅里叶变换，傅里叶变换提供了一个将信号从时域转变到频域的方法。之所以要有信号的频域分析，是因为很多信号在时域不明显的特征可以在频域下得到很好的展现，可以更加容易的进行分析和处理。
FFT
Matlab提供的傅里叶变换的函数是FFT，中文名叫做快速傅里叶变换。快速傅里叶变换的提出是伟大的，使得处理器处理数字信号的能力大大提升，也使我们生活向数字化迈了一大步。
接下来就谈谈如何使用这个函数。
fft使用很简单，但是一般信号都有x和y两个向量，而fft只会处理y向量，所以想让频域分析变得有意义，那么就需要用户自己处理x向量
一个简单的例子
从一个简单正弦信号开始吧，正弦信号定义为：
y(t)=2sin⁡(2πf0t)" role="presentation">y(t)=2sin(2πf0t)y(t)=2sin(2πf0t)
我们现在通过以下代码在Matlab中画出这个正弦曲线
这就是我们得到的：
当我们对这条曲线fft时，我们希望在频域得到以下频谱（基于傅里叶变换理论，我们希望看见一个幅值为1的峰值在-4Hz处，另一个在+4Hz处）
使用FFT命令
我们知道目标是什么了，那么现在使用Matlab的内建的FFT函数来重新生成频谱
1
2
3
4
5
6
7
8
9
10
11
12
13
%plot the frequency spectrum using the MATLAB fft command
matlabFFT = figure; %create a new figure
YfreqDomain = fft(y); %take the fft of our sin wave, y(t)
stem(abs(YfreqDomain)); %use abs command to get the magnitude
%similary, we would use angle command to get the phase plot!
%we'll discuss phase in another post though!
xlabel('Sample Number')
ylabel('Amplitude')
title('Using the Matlab fft command')
grid
axis([0,100,0,120])
效果如下：
但是注意一下，这并不是我们真正想要的，有一些信息是缺失的
x轴本来应该给我们提供频率信息，但是你能读出频率吗？
幅度都是100
没有让频谱中心为0
为FFT定义一个函数来获取双边频谱
以下代码可以简化获取双边频谱的过程，复制并保存到你的.m文件中
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
function [X,freq]=centeredFFT(x,Fs)
%this is a custom function that helps in plotting the two-sided spectrum
%x is the signal that is to be transformed
%Fs is the sampling rate
N=length(x);
%this part of the code generates that frequency axis
if mod(N,2)==0
k=-N/2:N/2-1; % N even
else
k=-(N-1)/2:(N-1)/2; % N odd
end
T=N/Fs;
freq=k/T; %the frequency axis
%takes the fft of the signal, and adjusts the amplitude accordingly
X=fft(x)/N; % normalize the data
X=fftshift(X); %shifts the fft data so that it is centered
这个函数输出正确的频域范围和变换后的信号，它需要输入需要变换的信号和采样率。
接下来使用前文的正弦信号做一个简单的示例，注意你的示例.m文件要和centeredFFT.m文件在一个目录下
1
2
3
4
5
6
7
8
9
10
[YfreqDomain,frequencyRange] = centeredFFT(y,Fs);
centeredFFT = figure;
%remember to take the abs of YfreqDomain to get the magnitude!
stem(frequencyRange,abs(YfreqDomain));
xlabel('Freq (Hz)')
ylabel('Amplitude')
title('Using the centeredFFT function')
grid
axis([-6,6,0,1.5])
效果如下：
这张图就满足了我们的需求,我们得到了在+4和-4处的峰值，而且幅值为1.
为FFT定义一个函数来获取右边频谱
从上图可以看出，FFT变换得到的频谱是左右对称的，因此，我们只需要其中一边就能获得信号的所有信息，我们一般保留正频率一侧。
以下的函数对上面的自定义函数做了一些修改，让它可以帮助我们只画出信号的正频率一侧
1
2
3
4
5
6
7
8
9
10
11
12
13
function [X,freq]=positiveFFT(x,Fs)
N=length(x); %get the number of points
k=0:N-1; %create a vector from 0 to N-1
T=N/Fs; %get the frequency interval
freq=k/T; %create the frequency range
X=fft(x)/N; % normalize the data
%only want the first half of the FFT, since it is redundant
cutOff = ceil(N/2);
%take only the first half of the spectrum
X = X(1:cutOff);
freq = freq(1:cutOff);
和前面一样，使用正弦信号做一个示例，下面是示例代码
1
2
3
4
5
6
7
8
9
[YfreqDomain,frequencyRange] = positiveFFT(y,Fs);
positiveFFT = figure;
stem(frequencyRange,abs(YfreqDomain));
set(positiveFFT,'Position',[500,500,500,300])
xlabel('Freq (Hz)')
ylabel('Amplitude')
title('Using the positiveFFT function')
grid
axis([0,20,0,1.5])
效果如下：
本文内容主要来自Matlab官方FFT教程，但是我找不到原链接了，就把我电脑上的文件上传到了百度云盘，公众号后台回复
下载|FFT_tutorial
获取下载链接
本文代码以上传github，可自行下载测试，链接：
https://github.com/greedyhao/DSP/tree/master/matlab_tutorial
