你的浏览器版本过低,可能导致网站不能正常访问!
为了你能正常使用网站功能,请使用这些浏览器。

【安富莱——DSP教程】第27章 FFT的Matlab实现

[复制链接]
baiyongbin2009 发布时间:2015-4-11 10:51
本帖最后由 baiyongbin2009 于 2015-4-11 10:58 编辑
2 r( R1 r, A7 ]4 i2 v
4 i+ f, @9 `/ @/ g3 ^' }特别说明:完整45期数字信号处理教程,原创高性能示波器代码全开源地址:链接
, P5 v: D2 v1 \! R; P" z9 W) G$ L+ X
第27章 FFT的Matlab实现

& u; r. d7 P$ T9 o8 O8 ]
    本章主要讲解fft,ifft和fftshift在matlab上的实现。
    27.1 FFT函数
    27.2 IFFT函数
    27.3 FFTSHIFT函数
    27.4 总结

3 n8 B% Q3 _, v5 M- q# O8 c, m% H27.1 FFT函数: h  y1 G. B7 i3 \: t
27.1.1 语法
Y = fft(x)
Y = fft(X,n)
Y = fft(X,[],dim)
Y = fft(X,n,dim)

- n2 l5 A9 M) \3 u
27.1.2 定义
Y = fft(x) 和 y = ifft(X)分别用于实现正变换和逆变换,公式描述如下:
27.1.png
' _8 `$ R6 F/ l- ~" [- G
27.1.3 描述
Y = fft(x)
    此函数用于返回向量x的离散傅立叶变换(DFT),计算时使用快速傅里叶算法(fast Fourier transform (FFT))。
    如果输入X是一个矩阵,Y = fft(X)返回该矩阵的每列的傅里叶变换。
    如果输入X是一个多维数组,实现第一个尺寸不为1的维度的FFT变化,注意这里第一个尺寸不为1是指一个矩阵的第一个尺寸不为1的维。
    比如一个矩阵是2*1,那么第一个尺寸不为1的维就是行(尺寸为2)
    X是 1*2*3表示第一个尺寸不为1的维就是列(尺寸为2)
    X为维数5*6*2的话,第一个尺寸不为1的维就是行(尺寸为5))
Y = fft(X,n)
    此函数用于返回n点的DFT。fft(n)和fft(X,n)是等同的,其中n是向量X中第一个尺寸不为1的维度。如果X的长度小于n,则X的长度通过填充零达到长度为n。如果X的长度大于n,序列X被截断。当X是一个矩阵,各列的长度都以相同的方式进行调整。
Y = fft(X,[],dim)
Y = fft(X,n,dim)
    上面两个函数用于实现指定维度的FFT运算。

2 l/ b! |4 w* H5 |
27.1.4 FFT实例一:幅频响应
    傅里叶变换的一个常见用途就是查找埋藏在噪声信号中的实际信号的频率成分。下面我们考虑一个这样的例子:
    采样率是1000Hz ,信号由如下三个波形组成。
      (1)50Hz的正弦波、振幅0,7。
      (2)70Hz正弦波、振幅1。
      (3)均值为0的随机噪声。
实际运行代码如下:
Fs = 1000;          %采样率
T = 1/Fs;           %采样时间单位
L = 1000;           %信号长度
t = (0-1)*T;        %时间序列
' f+ w7 Q* v6 O2 K9 R8 Y
x = 0.7*sin(2*pi*50*t) + sin(2*pi*120*t);  %原始信号
y = x + 2*randn(size(t));     %原始信号叠加了噪声后
plot(Fs*t(1:50),y(1:50));      %绘制波形
title('原始信号+零均值随机噪声 ');
xlabel('时间单位:ms');
运行Matlab后,显示波形如下:
27.2.png
通过上面的截图,我们是很难发现波形中的频率成分,下面我们通过FFT变换,从频域观察就很方便了,Matlab运行代码如下:
Fs = 1000;         %采样率
T = 1/Fs;           %采样时间单位
L = 1000;          %信号长度
t = (0-1)*T;       %时间序列
6 V! j8 ^* w  \& h. E1 J
x = 0.7*sin(2*pi*50*t) + sin(2*pi*120*t);  %原始信号
y = x + 2*randn(size(t));    %原始信号叠加了噪声后
& y" f3 y& c/ i% C- W& N2 ~1 o1 g! [3 u
NFFT = 2^nextpow2(L);  %求得最接近采样点的2^n,由于上面是1000点,那么最近的就是1024点。
Y=fft(y,NFFT)/L;          % 进行FFT变换,除以总的采样点数,方便观察实际值。                     
f = Fs/2*linspace(0,1,NFFT/2+1); %频率轴,这里只显示Fs/2部分,另一半是对称的。
9 @% n/ {6 p# j0 L8 I, [
plot(f,2*abs(Y(1:NFFT/2+1))) %绘制波形
title('幅频相应');
xlabel('频率');
ylabel('幅度');
27.3.png
从上面的幅频相应,我们可以看出,两个正弦波的频谱并不是准确的0.5和1,而是比较接近,这个就是咱们在上节教程中所示的频谱泄露以及噪声的干扰。

: N( G* ^3 h2 H# r+ v5 _7 D; b% Y
27.1.5 FFT实例二:相频响应
    这里我们以采样率=256Hz采样信号1.5*sin(2*pi*50*t+pi/3),并求出其幅频和相频响应,Matlab上面运行的代码如下:
Fs = 256;              % 采样率
N  = 256;             % 采样点数
n  = 0:N-1;            % 采样序列
t  = 0:1/Fs:1-1/Fs;     % 时间序列
f = n * Fs / N;          %真实的频率

) _% H2 c: D2 d
x = 1.5*sin(2*pi*50*t+pi/3) ;  %原始信号
y = fft(x, N);    %对原始信号做FFT变换
Mag = abs(y);    %求FFT转换结果的模值
subplot(2,1,1);
plot(f, Mag);       %绘制幅频相应曲线
title('幅频相应');
xlabel('频率/Hz');
ylabel('幅度');
' i. O. Z8 ]+ a  ~6 _
subplot(2,1,2);
plot(f,  angle(y)*180/pi);   %绘制相频响应曲线,注意这将弧度转换成了角度
title('相频响应');
xlabel('频率/Hz');
ylabel('幅度');
运行后求出的幅频相应和相频响应结果如下:
27.4.png

# q, {6 a% D9 I3 G" n" Q
收藏 评论14 发布时间:2015-4-11 10:51

举报

14个回答
baiyongbin2009 回答时间:2015-4-11 10:55:06
27.2 IFFT函数' A1 C3 H( o* G; ?' D
27.2.1 语法
y = ifft(X)
y = ifft(X,n)
y = ifft(X,[],dim)
y = ifft(X,n,dim)
y = ifft(..., 'symmetric')
y = ifft(..., 'nonsymmetric')

+ v1 L; d: G: X
27.2.2 描述
y = ifft(X)
    此函数用于返回向量X的离散傅立叶变换(DFT)逆变换结果,计算时使用快速傅里叶算法(fast Fourier transform (FFT))。
y = ifft(X,n)
    此函数用于返回n点的IDFT。
y = ifft(X,[],dim)
y = ifft(X,n,dim)
    上面两个函数用于实现指定维度的IFFT运算。
5 {5 Y6 |8 g. I7 D1 }+ A! P& h$ A
27.2.3 IFFT实例
    下面我们对信号:0.7*sin(2*pi*50*t) + sin(2*pi*120*t)求FFT和IFFT,并绘制原始信号和转换后的信号。Matlab上运行的代码如下:
Fs = 1000;             %采样率
T = 1/Fs;               % 采样时间
L = 1024;               % 信号长度
t = (0-1)*T;           % 时间序列
y = 0.7*sin(2*pi*50*t) + sin(2*pi*120*t);  %50Hz正弦波和120Hz的正弦波的叠加
subplot(2,1,1);
plot(Fs*t(1:50), y(1:50));   %绘制原始信号
title('原始信号');
Y = fft(y, L);               %分别调用正变换和逆变换                 
Z = ifft(Y);
subplot(2,1,2);
plot(Fs*t(1:50), Z(1:50));  %绘制逆变换后的波形
title('FFT和IFFT转换后的信号');
运行后求出的结果如下:
27.5.png
通过上面的运行结果可以看出,转换后的波形与原始的波形基本是一样的。

) C2 M0 X3 H% u7 p0 D
baiyongbin2009 回答时间:2015-4-11 10:57:28
27.3 FFTSHIFT函数
    fftshift的作用正是让正半轴部分和负半轴部分的图像分别关于各自的中心对称。因为直接用fft得出的数据与频率不是对应的,fftshift可以纠正过来
    以下是Matlab的帮助文件中对fftshift的说明:
    Y = fftshift(X) rearranges the outputs of fft, fft2, and fftn by moving the zero-frequency component to the center of the array. It is useful for visualizing a Fourier transform with the zero-frequency component in the middle of the spectrum. For vectors, fftshift(X) swaps the left and right halves of X。
    下面我们在Matlab上面实现一个如下的代码来说明fftshift的使用:
Fs = 256;              % 采样率
N  = 256;             % 采样点数
n  = 0:N-1;            % 采样序列
t  = 0:1/Fs:1-1/Fs;     % 时间序列
f = (-N/2:N/2-1) * Fs / N;   %真实的频率
x = 1.5*sin(2*pi*50*t+pi/3) ;  %原始信号
y = fft(x, N);      %对原始信号做FFT变换
Mag = abs(y);    %求FFT转换结果的模值
subplot(2,1,1);
plot(f, Mag);       %绘制幅频相应曲线
title('fft幅频相应');
xlabel('频率/Hz');
ylabel('幅度');
z = fftshift(y);    %对FFT转换后的结果做偏移。
subplot(2,1,2);
plot(f, z);        %绘制幅频相应曲线
title('fftshift幅频相应');
xlabel('频率/Hz');
ylabel('幅度');
Matlab的运行结果如下:
27.6.png
通过上面的运行结果我们可以看到,经过fftshift的调节后,正弦波的中心频率正好对应在了相应的50Hz频率点。使用fftshift还有很多其它的好处,有兴趣的可以查找相关的资料进行了解。
+ L' Y' q# \, m/ C
baiyongbin2009 回答时间:2015-4-11 10:58:16
27.4 总结
    本章节主要讲解了fft,iff和fftshift的基本用法,如果要深入了解,一定要多练习,多查资料和翻阅相关书籍。

8 U4 X" k8 i6 n1 v
kqh1120 回答时间:2015-4-11 11:46:15
学习了啊 smile.gif
eurphan 回答时间:2015-4-11 18:36:54
好东西   
zhangdaijin 回答时间:2015-4-11 23:28:42
学习学习
拼命三郎 回答时间:2015-4-12 10:50:55
学习学习,高大上的东西
拼命三郎 回答时间:2015-4-12 10:51:34
FFT好难啊
大器所成 回答时间:2015-4-12 13:12:25
有机会一定好好学学这FFT
wamcncn 回答时间:2015-4-12 23:00:51
标记下,学习
木木鱼 回答时间:2015-4-13 14:55:45
看后懂了很多
wyxy163@126.com 回答时间:2015-4-13 19:21:43
提示: 作者被禁止或删除 内容自动屏蔽
caiyuhui74748 回答时间:2015-7-29 09:42:44
好东西顶起来,学起来
烟花绽放 回答时间:2015-7-30 08:53:17
学习了,

所属标签

关于
我们是谁
投资者关系
意法半导体可持续发展举措
创新与技术
意法半导体官网
联系我们
联系ST分支机构
寻找销售人员和分销渠道
社区
媒体中心
活动与培训
隐私策略
隐私策略
Cookies管理
行使您的权利
官方最新发布
STM32N6 AI生态系统
STM32MCU,MPU高性能GUI
ST ACEPACK电源模块
意法半导体生物传感器
STM32Cube扩展软件包
关注我们
st-img 微信公众号
st-img 手机版