65.9K
CodeProject 正在变化。 阅读更多。
Home

使用 Matlab 进行简单的 FFT 和滤波教程

Grasshopper.iics
starIconstarIconstarIconstarIcon
emptyStarIcon
starIcon

4.27/5 (4投票s)

2012年1月20日

CPOL

4分钟阅读

viewsIcon

151532

理解 FFT 从未如此简单！

引言

让我们理解 FFT。它是快速傅里叶变换，是一种快速有效地计算 DFT 或离散傅里叶变换的算法。看到标题出现的第一个问题是，在 2012 年，当该算法已经有大约 50 年历史时，Code Project 的新文章部分中关于 FFT 的教程在做什么。 好问题。 顺便说一下，您可以参考以下链接，以复杂的方式理解相同的内容。

http://en.wikipedia.org/wiki/Discrete_Fourier_transform

http://en.wikipedia.org/wiki/Fast_Fourier_transform

背景

我认为我必须编写本教程的原因之一是，我看到许多工程师发现很难理解什么是频域，什么是滤波。 我们如何以简单的方式过滤信号。 什么是 FFT 的对称特征，更重要的是理解关于 FFT 的核心事实。 因此，如果您是一位编程高手，并且是一位编写了自己的 DSP 逻辑和算法的杀手级程序员，那么本文可能不会给您带来任何新的或有趣的东西。 但是，如果您刚刚开始接触信号处理，您可能会觉得它很有帮助（希望如此）。

Using the Code

首先，让我们构造一个幅度为 5 的 50Hz 简单正弦信号。 我考虑的采样频率是消息频率的 100 倍（它至少应该是奈奎斯特速率的两倍），这意味着我将从实际的模拟正弦信号中收集 1000 个样本。 因此，如果频率为 50，则意味着您将在 1/50 或 0.02 秒处获得一个样本。 如果您想要 10 个样本，则必须收集 0.2 秒的样本。 您将在每 1/5000 秒采样一次信号。

f=50;
A=5;
Fs=f*100;
Ts=1/Fs;
t=0:Ts:10/f;
x=A*sin(2*pi*f*t);
plot(x),grid on

1.png

图 1：频率为 50Hz 的正弦信号。 变量 x 的图。

所以这张图显示了 50 Hz 的 10 个正弦波周期。 您可以将绘图更改为 stem 类型以查看样本。

那么让我们进行 FFT。

F=fft(x);

您一定已经学习过 N 点 FFT。 如果您不指定 N，则默认情况下 N 是消息信号的长度。 在本例中，它是 1001。

非常重要的事情：FFT 将您的采样频率划分为 N 个相等的部分，并返回每个频率级别的信号强度。 这意味着您将频率从 0 到 5000 划分为 1001 个相等的部分。 因此，fft 中的第一个点是 5Hz，下一个表示 10 Hz，依此类推。 由于您的消息的实际频率为 50 Hz，因此您必须在该级别获得最高峰值。 那么让我们绘制 FFT。

plot(real(F)),grid on

2.png

图 2：信号 x 的 FFT 输出实部的图，并标出了最高峰值。 变量 F 实部的图。

你看到了什么？ 最高值位于第 11 个位置，这不过是第 10 个样本，因为 matlab 从 1 开始。

10*5?=50Hz

这太简单了。 但是等等，为什么我们也会得到虚部？ 绝对没有什么叫做虚数值，它们反映了过去的值。

您的信号是正弦信号，所以您也会得到负周期，对吧。 检查这个。

plot(real(imag(F))),grid on

3.png

图 3：信号 x 的 FFT 输出虚部的图，并标出了最高峰值。 变量 F 虚部的图。

你看到了什么？ :)

是的，同样的事情，也有一个负幅度峰值。 我希望你明白了。

现在让我们了解如何使用 FFT 进行一些滤波。 我要做的是将几个具有不同频率和相同样本数的信号与您的信号 x 混合。

x1=A*sin(2*pi*(f+50)*t);
x2=A*sin(2*pi*(f+250)*t);
x=x+x1+x2;

4.png

图 4：包含 50Hz、100Hz、300Hz 频率分量的混合信号 x 的图。 变量 x 的图。

我们的信号现在包含三个分量：50Hz、100Hz 和 300Hz。

它的 fft 图就是证明。 三个频率在三个级别，正好在 50、100 和 300Hz。

5.png

图 5：混合信号 x 的 FFT 输出实部的图。

我想恢复我的原始信号。 因此，如果我将样本保留在第 11 个样本处，而将其余样本保留，然后进行 ifft，我必须取回我的信号，对吗？ 但是必须保留第 10 个和第 11 个样本的值。 这称为加窗。

6.png

图 6：在混合信号 x 的 FFT 图中，将所需频带标记为矩形窗口。

请参阅我在混合信号的 FFT 响应上放置的红色窗口。

F2=zeros(length(F),1);
F2(10:11)=F(10:11);
xr=ifft(F2);
plot(real(xr)),grid on

7.png

图 7：使用矩形窗口从图 4 所示的混合信号中恢复频率分量为 50Hz 的原始信号。

你看到了什么？ 这是您的原始信号。 但是幅度较小。 这是因为您选择了矩形窗口，并且损耗是由于负滤波器增益造成的。 您现在可以阅读一些数字滤波器并应用它。

这是完整的代码，可能不需要，但我仍然发布它，以防您决定花一些时间并使用它。

 clear all;
clc;
close all
f=50;
 A=5;
Fs=f*100;
Ts=1/Fs;
t=0:Ts:10/f;
x=A*sin(2*pi*f*t);
x1=A*sin(2*pi*(f+50)*t);
x2=A*sin(2*pi*(f+250)*t);
x=x+x1+x2;
plot(x)
F=fft(x);
figure
N=Fs/length(F);
baxis=(1:N:N*(length(x)/2-1));
plot(baxis,real(F(1:length(F)/2)))

关注点

您可以从此来源学习 Matlab 基础知识 <这里>

要了解有关任何 Matlab 命令的详细信息，您可以简单地单击编辑器中的该命令并按 F1。 Matlab 帮助文件解释了诸如 fft、sin 等命令的用法和其他详细信息。

在所有复杂的数学背后，都有一个简单的逻辑。 我真的很喜欢那个简单的逻辑。

© . All rights reserved.