实验笔记【3】| 利用matlab进行fft分析

FiveMay

已于 2024-06-30 12:24:59 修改

阅读量1.5k

点赞数 4

文章标签：笔记

于 2024-06-30 12:21:15 首次发布

本文链接：https://blog.csdn.net/weixin_45884541/article/details/139700880

版权

主体内容来源：MATLAB中利用FFT进行信号频谱分析：原理、技术与实践-百度开发者中心 (baidu.com)

在信号处理中，频谱分析是一种关键技术，用于理解信号的频率成分。快速傅里叶变换（FFT）是实现频谱分析的常用方法，特别是在处理实时信号时。在MATLAB中，我们可以利用FFT函数来进行频谱分析。下面，我们将分步骤解释如何使用FFT在MATLAB中进行信号频谱分析。

一、FFT的基本原理

FFT是一种高效的算法，用于计算离散傅里叶变换（DFT）和其逆变换。通过将信号分解成多个频率分量，我们可以了解信号的频率特性。在FFT中，我们首先将信号分成多个更短的信号段（通常称为帧），然后对每个段执行DFT。

二、选择合适的FFT大小

FFT的大小（N）通常等于帧的长度。对于给定的采样率，如果我们要分析的最高频率是f{max}，则合适的FFT大小应该是N=2^k，其中k是满足f{max} * 2π >= N/2的正整数。这是因为频率范围被限制在0到f_{max}之间，而N/2代表了最大的频率分量。

三、处理频谱泄漏和窗函数

频谱泄漏是指在DFT过程中产生的伪频率分量。这些是由窗口函数的非理想特性引起的。为了减小频谱泄漏，可以使用不同类型的窗函数（如汉宁窗、哈明窗等）。在MATLAB中，我们可以使用’hamming’或’hanning’等选项来选择不同的窗函数。

四、代码示例

下面是一个简单的MATLAB代码示例，展示如何使用FFT进行频谱分析：

% 创建一个包含两个频率分量的信号
Fs = 1000; % 采样率
t = 0:1/Fs:1-1/Fs; % 时间向量
x = cos(2*pi*50*t) + cos(2*pi*120*t); % 合成信号
% 使用FFT进行频谱分析
N = length(x); % 获取信号长度
X = fft(x, N)/N; % 计算FFT并归一化
f = Fs*(0:(N/2))/N; % 频率向量
% 绘制频谱图
plot(f, abs(X(1:N/2+1))) % 只绘制一半的频谱，因为X是对称的
title('频谱分析')
xlabel('频率 (Hz)')
ylabel('幅度')

这段代码首先创建了一个包含两个频率分量的合成信号。然后，它使用FFT计算信号的频谱，并绘制出频谱图。注意，由于FFT的结果是对称的，我们通常只绘制一半的频谱。

五、结论

通过以上步骤，我们可以看到在MATLAB中使用FFT进行信号频谱分析的完整过程。理解FFT的基本原理以及如何选择合适的FFT大小是非常重要的。同时，处理频谱泄漏和窗函数也是关键的考虑因素。通过实际代码示例，我们可以更好地理解如何将这些概念应用于实际问题中。希望本文能帮助您更好地理解和应用频谱分析技术。

六、补充内容（频谱泄露VS窗函数）

cr：深入浅出的理解频谱泄露-CSDN博客

1）什么是频谱泄露？

　　频谱泄露与 傅里叶变换尤其是离散时间傅里叶变换有关，对于频谱泄露，通常的解释是这样的：

　　信号为无限长序列，运算需要截取其中一部分（截断），于是需要加窗函数，加了窗函数相当于时域相乘，于是相当于频域卷积（两个信号频域卷积后最高频率是两个频率相加），于是频谱中除了本来该有的主瓣之外，还会出现本不该有的旁瓣，这就是频谱泄露！为了减弱频谱泄露，可以采用加权的窗函数，加权的窗函数包括平顶窗、汉宁窗、高斯窗等等。而未加权的矩形窗泄露最为严重。

为了说明频谱泄露，一下子引入了时域、频域、窗函数、卷积、主瓣、旁瓣等等抽象的概念。频谱泄露有这么复杂吗？频谱泄露到底是什么意思？

　　一句话，频谱泄露就是分析结果中，出现了本来没有的频率分量。比如说，50Hz的纯正弦波，本来只有一种频率分量，分析结果却包含了与50Hz频率相近的其它频率分量。

　　更简单的描述是：分析结果与实际不一致！

2）为何会出现频谱泄露？

　　我们把无限长序列分为两种情况：

1、无限长序列为非周期信号

　　非周期的无限长序列，任意截取一段有限长的序列，都不能代表实际信号，分析结果当然与实际信号不一致！

2、无限长序列为周期信号

　　周期的无限长序列，假设截取的是正好一个或整数个信号周期的序列，这个有限长序列就可以代表原无限长序列，如果分析的方法得当的话，分析结果应该与实际信号一致！

　　第一种情况，我们可以这样理解：你分析的内容根本就不能代表实际内容，结果当然不一致，更准确的说法，结果是错误的，造成错误结果的原因是分析方法是错误的！　

第二种情况，我们作了两个假设，第二个假设是成立的，该假设基于伟大的傅里叶作出的论断！那么，如果第一个假设也成立，是不是就不会发生频谱泄露呢？

　　答案是肯定的！

　　从无限长序列中截取一个或整数个周期，我们称为整周期截断，反之，称为 非整周期截断。

　　 整周期截断，不会造成频谱泄露！

　　非整周期截断，必然造成频谱泄露！

　　换言之：

　　 整周期截断是不发生频谱泄露的充分且必要条件！

　　或：

　　 非整周期截断是发生频谱泄露的充分且必要条件！

　　为什么非整周期截断就会发生频谱泄露呢？

原因实际上与第一种情况是类似的。如图1所示，截取分析的有限长序列，傅里叶变换仍然将信号当成无限长序列，对于有限长序列，又是如何当成无限长呢？

　　这里采用了一种被称为周期延拓的技术，所谓周期延拓，就是把截取的有限长序列当成是无限长序列的一个周期，然后不断的复制，得到一个新的无限长序列。

如图2所示，从图1所示无限长序列中截取的有限长序列，经过周期延拓后，得到一个新的无限长序列，显然，这个新的序列与原序列是不一样的！

　　图2的信号与图1的信号不同，分析得到的频谱自然也不同！不同之处在于，图1是单一频率信号，只有一根谱线，而图2中，除了图1信号包含的这根谱线（不妨称为主谱线）外，出现了其它频率的谱线，通常，这些谱线要比主谱线短很多，如果把这些原信号不包含的谱线理解为是主谱线泄露出来的，那么，这种现象就被称为频谱泄露！

　　采用合适的窗函数（常见的窗函数有汉宁窗、三角窗、海明窗和高斯窗等等）可以一定程度上抑制频谱泄露。

　　窗函数的概念，非常抽象，然而，窗函数的作用，是非常有限的，我们可以这样理解：

　　如图2中的信号，由于突然截断造成周期延拓时两个周期相邻处出现了信号突变，这种突变，代表的是信号包含了高次谐波。加上合适的歘窗函数，可以把这个突变变得圆滑一些，从而抑制高次谐波。

　　但是，我们也可以这样想，假设图2的信号就是真实信号，那么，加上这样的窗函数反而得到了错误的结果！

　　因此，避免频谱泄露的根本还是要从源头出发，尽可能做到准确的整周期截断，这种情况下，窗函数可以选择最简单的矩形窗。

3）正确处理频谱泄露

　　造成频谱泄露的原因在于傅里叶变换的输入信号不能准确的、完整的代表被分析信号，输出产生的一种误差，这种误差可以通过加合适的窗函数或延长时间窗得以改善，当输入信号的不完整性达到一定程度，输出是一种错误的结果。

　　对于周期信号，整周期截断是不发生频谱泄露的充分且必要条件，抑制频谱泄露应该从源头抓起，尽可能进行整周期截断。

七、波特图

1、伯德图的定义

伯德图由对数幅频特性和对数相频特性两条曲线构成。

伯德图是线性非时变系统的传递函数对频率的半对数坐标图，其横轴频率以对数尺度（log scale）表示，纵坐标幅值或相角采用线性分度，利用伯德图可以看出系统的频率响应。

伯德图一般是由二张图组合而成，伯德图由两张图组成：①G(jω)的幅值(以分贝，dB表示)-频率(以对数标度)对数坐标图，其上画有对数幅频曲线；②G(jω)的相角-频率(以对数标度)对数坐标图，其上画有相频曲线。

对数幅值的标准表达式为20 lg|G(jω)|，单位是分贝，相角的单位是度。由于增益用对数来表示（log(ab)=log(a)+log(b)），因此一传递函数乘以一常数，在伯德增益图只需将图形的纵向移动即可，二传递函数的相乘，在波德幅频图就变成图形的相加。幅频图纵轴0分贝以下具有正增益裕度、属稳定区，反之属不稳定区。

配合波德相频图可以估算一信号进入系统后，输出信号及原始信号的比例关系及相位。例如一个Asin(ωt) 的信号进入系统后振幅变原来的k倍，相位落后原信号Φ，则其输出信号则为(Ak)sin(ωt−Φ)，其中的k和Φ都是频率的函数。相频图纵轴-180度以上具有正相位裕度、属稳定区，反之属不稳定区。

2、伯德图的画法

画伯德图时，分三个频段进行，先画幅频特性，顺序是中频段、低频段和高频段。将三个频段的频率特性(或称频率响应)合起来就是全频段的幅频特性，然后再根据幅频特性画出相应的相频特性来。

作伯德图时，首先写出频率特性，然后按常数因子K、积分和微分因子(jω)、一阶因子(1+jωT)和二阶因子[1+2ζ(jω/ωn)+(jω)/ω]1这样四种基本因子分别画出伯德图，再总加而成。

下图是常见简易传递函数的伯德图趋势：

在matlab中通常绘制语句为semilogx(f,20*log10(P1))

p.s: matlab读取数据的函数？

csvread() 函数-----读取逗号分隔的纯数值csv文件

csv文件就是comma-separated value (CSV) file。数据使以逗号相隔的形式保存在.csv文件中。
2019最新版的官方文档不在推荐使用csvread读取csv文档，而是推荐使用readmatrix但是目前还是兼容。
三种方式读取：
(1) M = csvread(filename), 文件的内容只能是数值。
(2) M = csvread(filename,R1,C1)，指定从R1+1行与C1+1列开始读其后的所有内容。
(3) M = csvread(filename,R1,C1,[R1 C1 R2 C2])，通过指定左上角开始的行列和右下角的行列读取的范围。

注：① 方法一读取整个文件，必须保证该文件只有数值内容，一般用不到这种方法，毕竟从示波器等导出的文件中，一定会有描述内容，即各种单位等头文件，都是英文。

②R和C从零开始，因此R=0和C=0指定文件中的第一个值。即csv文件中的第一行为该函数的第0行，csv文件中的第一列为该函数的第0列；即从R行开始，C列开始，一直到结束的这个范围内的数字，全部存入M是一个矩阵，这种方式比较常用。

③（R1，C1）是要读取的数据的左上角，（R2，C2）是右下角。也可以使用电子表格表示法指定RNG，如RNG=’A1..B7’；从R1行C1列到R2行C2列之间读取数据。

④csvread()只能读取完全由逗号分隔的数字，矩阵中含有用来分隔数字的分号字符，会将这个矩阵当作字符串来看待；