关于这个问题,在很早之前就分享过,也通过了解实现了算法,当时看的明白,想的明白,突然要用的时候,又开始疑问,不免有些纠结,与其每次使用的时候都查,浪费时间,还不如,一次搞定。

真心没把哪门没学好的课程,归结到老师,但fft这事,还真得跟大学老师讨个说法,哈哈。

matlab FFT 横坐标问题:前人关于FFT横坐标的详细阐述

我们知道Fourier分析是信号处理里很重要的技术,matlab提供了强大的信号处理能力,但是有一些细节部分需要我们注意。

记信号f(t)的起始时间为t_start, 终止时间为t_end, 采样周期为t_s, 可以计算信号的持续时间Duration为 t_end – t_start, 信号离散化造成的采样点数 N = Duration/t_s + 1;

根据Fourier分析的相关结论,我们知道时域的采样将会造成频域的周期化,该周期为采样频率f_s(著名的香农采样定理基于此).

于是, 经过matlab的fft函数处理后,得到数据的横坐标为0:f_s/(N-1):f_s。相关代码如下所示:

%matlab fft 测试代码

t_s = 0.01;

t_start = 0.5; t_end = 5;

t = t_start:t_s:t_end;

y = 0.5*sin(2*pi*15*t)+2*sin(2*pi*40*t);

y_f = fft(y);

subplot(3,1,1);

plot(t,y); title('original signal');

Duration = t_end - t_start;

Sampling_points = Duration/t_s + 1;

f_s = 1/t_s;

f_x = 0:f_s/(Sampling_points-1):f_s;

subplot(3,1,2);

plot(f_x,abs(y_f)); title('fft transform');

subplot(3,1,3);

plot(f_x-f_s/2,abs(fftshift(y_f))); title('shift fft transform');

也就是说,如果我们不使用fftshift,其变换后的横坐标为0:f_s/(N-1):f_s,如果使用fftshift命令,0频率分量将会移到坐标中心,这也正是matlab中帮助中心给出的意思:对fft的坐标进行了处理。实际上由于频谱的周期性,我们这样做是合理的,可以接受的。

请读者特别要注意横坐标的差别。另外,根据函数的特性,频谱应当只有在15Hz,40Hz出现峰值,但是fft变换后在60Hz,及85Hz处同样出现了峰值,应当可以从fft的计算过程中得到相应的解释。

事实上,如果我们用15Hz,60Hz来测试fft变换,也即是 y = 0.5*sin(2*pi*15*t)+2*sin(2*pi*60*t);图像如下所示,没有任何变化。

这种现象提醒我们,频率在f_s以内,即 0

实际上,这也就间接地证明了Nyquist采样定理的合理性:采样频率要高于截止频率的两倍,上面的处理中我们所使用的采样频率为100Hz,于是当截止频率超过50Hz时,就会出现混叠效应,特殊情况就如上图所示:完全一样。于是,这也就告诉我们若要正确的显示频谱,需要仔细地考量采样频率与截止频率的关系,若太小,则有可能出现混叠,若太大,则计算代价过高。

上面内容说明了fft横坐标与采样定理的问题,如果你没有这个耐心看完,想着这事啥呀,乱七八糟的。下面我总结了坐标问题如下:

1)如果你用fft(NFFT)计算振幅谱,且fft变换后的点数为NFFT的话,只用fft,而不用fftshift,那么:

你只需要显示fft变换结果的前NFFT/2个点,横坐标为fx=index*Fs/(nfft-1),index=0:nfft/2-1 ;

2)如果你用fft(NFFT)计算振幅谱,且fft变换后的点数为NFFT的话,用了fft,还用了fftshift,那么:

你需要显示fft变换结果的数据点,横坐标为fx = 0:Fs/(nfft-1):Fs;fx=fx-Fs/2;

原因:对实信号fft,产生了负频率,由于FFT变化是对称的,所以负频率对应到了Fs/2~Fs上,对实信号而言,该段频率没有任何意义,所以只显示到Fs/2。

matlab FFT 纵坐标问题:前人关于FFT纵坐标的详细阐述

一.调用方法

X=FFT(x);

X=FFT(x,N);

x=IFFT(X);

x=IFFT(X,N)

用MATLAB进行谱分析时注意:

(1)函数FFT返回值的数据结构具有对称性。

例:

N=8;

n=0:N-1;

xn=[4 3 2 6 7 8 9 0];

Xk=fft(xn)

Xk =

39.0000 -10.7782 + 6.2929i 0 - 5.0000i 4.7782 - 7.7071i 5.0000 4.7782 + 7.7071i 0 + 5.0000i -10.7782 - 6.2929i

Xk与xn的维数相同,共有8个元素。Xk的第一个数对应于直流分量,即频率值为0。

(2)做FFT分析时,幅值大小与FFT选择的点数有关,但不影响分析结果。在IFFT时已经做了处理。要得到真实的振幅值的大小,只要将得到的变换后结果乘以2除以N即可。

二.FFT应用举例

例1:x=0.5*sin(2*pi*15*t)+2*sin(2*pi*40*t)。采样频率fs=100Hz,分别绘制N=128、1024点幅频图。

clf;

fs=100;N=128; %采样频率和数据点数

n=0:N-1;t=n/fs; %时间序列

x=0.5*sin(2*pi*15*t)+2*sin(2*pi*40*t); %信号

y=fft(x,N); %对信号进行快速Fourier变换

mag=abs(y); %求得Fourier变换后的振幅

f=n*fs/N; %频率序列

subplot(2,2,1),plot(f,mag); %绘出随频率变化的振幅

xlabel('频率/Hz');

ylabel('振幅');title('N=128');grid on;

subplot(2,2,2),plot(f(1:N/2),mag(1:N/2)); %绘出Nyquist频率之前随频率变化的振幅

xlabel('频率/Hz');

ylabel('振幅');title('N=128');grid on;

%对信号采样数据为1024点的处理

fs=100;N=1024;n=0:N-1;t=n/fs;

x=0.5*sin(2*pi*15*t)+2*sin(2*pi*40*t); %信号

y=fft(x,N); %对信号进行快速Fourier变换

mag=abs(y); %求取Fourier变换的振幅

f=n*fs/N;

subplot(2,2,3),plot(f,mag); %绘出随频率变化的振幅

xlabel('频率/Hz');

ylabel('振幅');title('N=1024');grid on;

subplot(2,2,4)

plot(f(1:N/2),mag(1:N/2)); %绘出Nyquist频率之前随频率变化的振幅

xlabel('频率/Hz');

ylabel('振幅');title('N=1024');grid on;

运行结果:

fs=100Hz,Nyquist频率为fs/2=50Hz。整个频谱图是以Nyquist频率为对称轴的。并且可以明显识别出信号中含有两种频率成分:15Hz和40Hz。由此可以知道FFT变换数据的对称性。因此用FFT对信号做谱分析,只需考察0~Nyquist频率范围内的福频特性。若没有给出采样频率和采样间隔,则分析通常对归一化频率0~1进行。另外,振幅的大小与所用采样点数有关,采用128点和1024点的相同频率的振幅是有不同的表现值,但在同一幅图中,40Hz与15Hz振动幅值之比均为4:1,与真实振幅0.5:2是一致的。为了与真实振幅对应,需要将变换后结果乘以2除以N。

例2:x=0.5*sin(2*pi*15*t)+2*sin(2*pi*40*t),fs=100Hz,绘制:

(1)数据个数N=32,FFT所用的采样点数NFFT=32;

(2)N=32,NFFT=128;

(3)N=136,NFFT=128;

(4)N=136,NFFT=512。

clf;fs=100; %采样频率

Ndata=32; %数据长度

N=32; ?T的数据长度

n=0:Ndata-1;t=n/fs; %数据对应的时间序列

x=0.5*sin(2*pi*15*t)+2*sin(2*pi*40*t); %时间域信号

y=fft(x,N); %信号的Fourier变换

mag=abs(y); %求取振幅

f=(0:N-1)*fs/N; %真实频率

subplot(2,2,1),plot(f(1:N/2),mag(1:N/2)*2/N); %绘出Nyquist频率之前的振幅

xlabel('频率/Hz');ylabel('振幅');

title('Ndata=32 Nfft=32');grid on;

Ndata=32; %数据个数

N=128; ?T采用的数据长度

n=0:Ndata-1;t=n/fs; %时间序列

x=0.5*sin(2*pi*15*t)+2*sin(2*pi*40*t);

y=fft(x,N);

mag=abs(y);

f=(0:N-1)*fs/N; %真实频率

subplot(2,2,2),plot(f(1:N/2),mag(1:N/2)*2/N); %绘出Nyquist频率之前的振幅

xlabel('频率/Hz');ylabel('振幅');

title('Ndata=32 Nfft=128');grid on;

Ndata=136; %数据个数

N=128; ?T采用的数据个数

n=0:Ndata-1;t=n/fs; %时间序列

x=0.5*sin(2*pi*15*t)+2*sin(2*pi*40*t);

y=fft(x,N);

mag=abs(y);

f=(0:N-1)*fs/N; %真实频率

subplot(2,2,3),plot(f(1:N/2),mag(1:N/2)*2/N); %绘出Nyquist频率之前的振幅

xlabel('频率/Hz');ylabel('振幅');

title('Ndata=136 Nfft=128');grid on;

Ndata=136; %数据个数

N=512; ?T所用的数据个数

n=0:Ndata-1;t=n/fs; %时间序列

x=0.5*sin(2*pi*15*t)+2*sin(2*pi*40*t);

y=fft(x,N);

mag=abs(y);

f=(0:N-1)*fs/N; %真实频率

subplot(2,2,4),plot(f(1:N/2),mag(1:N/2)*2/N); %绘出Nyquist频率之前的振幅

xlabel('频率/Hz');ylabel('振幅');

title('Ndata=136 Nfft=512');grid on;

结论:

(1)当数据个数和FFT采用的数据个数均为32时,频率分辨率较低,但没有由于添零而导致的其他频率成分。

(2)由于在时间域内信号加零,致使振幅谱中出现很多其他成分,这是加零造成的。其振幅由于加了多个零而明显减小。

(3)FFT程序将数据截断,这时分辨率较高。

(4)也是在数据的末尾补零,但由于含有信号的数据个数足够多,FFT振幅谱也基本不受影响。

对信号进行频谱分析时,数据样本应有足够的长度,一般FFT程序中所用数据点数与原含有信号数据点数相同,这样的频谱图具有较高的质量,可减小因补零或截断而产生的影响。

例3:x=cos(2*pi*0.24*n)+cos(2*pi*0.26*n)

(1)数据点过少,几乎无法看出有关信号频谱的详细信息;

(2)中间的图是将x(n)补90个零,幅度频谱的数据相当密,称为高密度频谱图。但从图中很难看出信号的频谱成分。

(3)信号的有效数据很长,可以清楚地看出信号的频率成分,一个是0.24Hz,一个是0.26Hz,称为高分辨率频谱。

可见,采样数据过少,运用FFT变换不能分辨出其中的频率成分。添加零后可增加频谱中的数据个数,谱的密度增高了,但仍不能分辨其中的频率成分,即谱的分辨率没有提高。只有数据点数足够多时才能分辨其中的频率成分。

======================评论======================================

1)FFT的算法与奈奎斯特频率没有关系,奈奎斯特频率只是用来避免频率的混叠。

2)对称是因为fft本身算法产生的吧!奈奎斯特频率是系统最高频率*2得到的最低采样频率,和fft没有关系啊,如果低于奈奎斯特频率采样,就会发生频谱混叠,是把大于系统最高频率的部分对称加到小于的部分,发生混叠现象。

3)"频谱图是以Nyquist频率为对称轴的"这句话不对。对称轴是N/2

上面内容说明了fft纵坐标的问题,如果你没有这个耐心看完,想着这事啥呀,乱七八糟的。下面我总结了坐标问题如下:

1)如果你用fft(NFFT)计算振幅谱,且fft变换后的点数为NFFT的话,那么将得到的振幅谱abs(fft(x))的结果乘以2然后除以N。

2)无。

标注横纵坐标MATLAB,matlab FFT 横纵坐标相关推荐

  1. DFT的计算、FFT的基础代码、FFT的横纵坐标问题(matlab)

    FFT的定义 FFT:快速傅里叶变换,是DFT的快速算法. DFT(Discrete Fourier Transform):离散傅里叶变换.在DTFT之后,将傅里叶变换的结果也进行离散化,就是DFT. ...

  2. Matlab中fft作频谱横纵坐标

    关于这个问题,在很早之前就分享过,也通过了解实现了算法,当时看的明白,想的明白,突然要用的时候,又开始疑问,不免有些纠结,与其每次使用的时候都查,浪费时间,还不如,一次搞定. 真心没把哪门没学好的课程 ...

  3. Matlab plotyy画双纵坐标图实例

    转载自:http://blog.sina.com.cn/s/blog_49d955150100lxoe.html Matlab plotyy画双纵坐标图实例 x = 0:0.01:20; y1 = 2 ...

  4. matlab 3个纵坐标,[转载]Matlab plotyy画双纵坐标图实例

    Matlab plotyy画双纵坐标图实例 x = 0:0.01:20; y1 = 200*exp(-0.05*x).*sin(x); y2 = 0.8*exp(-0.5*x).*sin(10*x); ...

  5. 示波器数据用matlab进行fft,示波器CSV波形数据导入Matlab进行FFT分析.doc

    示波器CSV波形数据导入Matlab进行FFT分析 1,将CSV文件拖到workspace窗口,弹出的Import Wizard窗口中,点选"Next",新窗口中选第二项" ...

  6. matlab示波器导出csv数据,示波器CSV波形数据导入Matlab进行FFT分析

    示波器CSV波形数据导入Matlab进行FFT分析 1,将CSV文件拖到workspace窗口,弹出的Import Wizard窗口中,点选"Next",新窗口中选第二项" ...

  7. Matlab 计算 FFT 的方法及幅值问题

    欢迎转载,但请一定要给出原文链接,标注出处,支持原创! 谢谢~ https://blog.csdn.net/qq_29225913/article/details/105467006 目录 1.Mat ...

  8. MATLAB中FFT的使用方法(频谱分析)

    原文地址:MATLAB中FFT的使用方法(频谱分析)作者:飞鸿 说明:以下资源来源于<数字信号处理的MATLAB实现>万永革主编 一.调用方法 X=FFT(x): X=FFT(x,N): ...

  9. 离散傅里叶变换DFT与FFT,MATLAB的FFT函数使用(原创)——如何使用fft()绘制出真正的频谱图像

    以前一直对MATLAB中fft()函数的使用一直存在疑惑,为什么要加一 些参数,并且如何确定这些参数,也查了许多资料,但很多都感觉只是 表面一说根本没有讲清其本质.但随着学习的推进,慢慢有所领悟,所 ...

最新文章

  1. 世界首富贝佐斯将“退休”
  2. 实验六 html网页设计,网页设计.html · 谢泽华/面向对象与软件工程实验二:网页模仿 - Gitee.com...
  3. 为LUKS加密的磁盘/分区做增量备份
  4. php获取当前整点时间_8.PHP的日期和时间
  5. LUA: lua基础.
  6. [转]Oh My Zsh,安装,主题配置
  7. spring aop 会根据实际情况(有无接口)自动选择 两种 动态代理(jdk和cglib)之一...
  8. 375 Inscribed Circles and Isosceles Triangles 等腰三角形 内接圆 圆周率PI表示
  9. oracle将字符串转成数组_【算法打卡】上升下降字符串
  10. celery 可视化_Celery部署爬虫(三)
  11. html 字符串 放到webbrowser,delphi 直接将html字符串读入WebBrowser中
  12. 《财经》杂志:盛大新浪梦纪实(完全版)
  13. Java+spring 基于ssm的幼儿园管理系统程序#毕业设计
  14. InnoDB关键特性之插入缓冲
  15. 使用晨曦账本记录店铺收支
  16. centos7解压rar包
  17. 网线属于计算机网络的哪一层,网线的种类分哪几种?
  18. python提取excel表格数据
  19. 爬虫学习day1遇到问题汇总(带参数的访问百度,代理,金山翻译只能翻译一个固定单词
  20. 关于计算机组装ppt,ENBM_PPT_chap01_V1.3 计算机组装.ppt

热门文章

  1. Spark:Jieba对数据库里提取的记录进行中文分词
  2. SwAV:Unsupervised Learning of Visual Features by Contrasting Cluster Assignments
  3. cad 激活出错解决办法
  4. 笔记本深度学习训练散热实用指南
  5. 如果地球是方的,环球旅行该怎么规划?
  6. 服务器拒绝发送文件怎么办,服务器拒绝了您发送离线文件
  7. c语言两个指针相等,C语言之指针,便于理解
  8. 万维网(World Wide Web)
  9. Excel删除空白行(亲自实践)
  10. 【实战】使用 Web Animations API 实现一个精确计时的时钟