matlab时频分析之短时傅里叶变换 spectrogram_matlab短时傅里叶变换-程序员宅基地

技术标签: matlab  短时傅里叶变换  时频分析  STFT  信号处理  

matlab时频分析之短时傅里叶变换 spectrogram

短时傅里叶变换常用于缓慢时变信号的频谱分析,可以观察沿时间变化的频谱信号。

其优点如下图所示,弥补了频谱分析中不能观察时间的缺点,也弥补了时域分析不能获取频率的缺点。
在这里插入图片描述

1 STFT的基本原理

基本原理可以理解为对一段长信号,截取每一段时间的短信号做fft,将得到的频谱图时间沿时间轴排列,及可得到时频的云图。
在这里插入图片描述

2 matlab中实现

这里采用最基础的spectrogram函数作为分析函数。虽然在其它很多地方都能找到能实现更高级用法的函数,比如逆变换:
https://ww2.mathworks.cn/matlabcentral/fileexchange/45577-inverse-short-time-fourier-transform-istft-with-matlab
但是这里只是单纯的做可视化分析。

代码如下:

%短时傅里叶变换展示
fs=2^10;    %采样频率 fs=65536hz
dt=1/fs;    %时间精度
timestart=-4;
timeend=4;
t=(0:(timeend-timestart)/dt-1)*dt+timestart;
L=length(t);

%设置信号
z=sin(2*pi*5.*t).*(t<-2)+sin(2*pi*10.*t).*(t>=-2&t<0)+...
    sin(2*pi*20.*t).*(t>=0&t<2)+sin(2*pi*40.*t).*(t>=2);
z2=wextend(1,'sym',z,round(length(z)/2));%镜像延拓

wlen=512;%设置窗口长度。窗口越长时间分辨率越差,频率分辨率越好。
hop=1;%每次平移的步长,最小为1。越小图像时间精度越好,但计算量大。
z2=wkeep1(z2,L+1*wlen);%中间截断


%做短时傅里叶
h=hamming(wlen);%设置海明窗的窗长
f=1:0.5:60;%设置频率刻度

[tfr2,f,t2]=spectrogram(z2,h,wlen-hop,f,fs);
tfr2=tfr2*2/wlen*2;
figure
imagesc(t2+timestart-wlen/fs/2,f,abs(tfr2))

3 窗口长度与分辨率

短时傅里叶变换的分辨率和窗口长度wlen的选取,频率的刻度f、平移的步长hop都有关系。

spectrogram函数可以选取较密的频率刻度来增加频域的分辨率,减小平移步长增加时间的分辨率。

但是对于短时傅里叶变换,最重要的还是窗口长度的选取。当频域刻度和平移步长足够密时,增加的只是生成图像的大小,但是物理层面的分辨率却没有改变。改变物理层面分辨率的就是窗口长度。

下图为同一个信号,选取不同窗口长度的图像。其余量不变。
在这里插入图片描述
窗口长度大的话,比如2048,频率能够清晰的显示出来,但是时间上有些模糊。窗口长度比较小比如256,时间上的分界线能够清晰显示,但是频率的值却不能清晰读出。

这就是窗口的影响。可以想象,当窗口大到整个信号长度时,短时傅里叶变换就退化为傅里叶变换,没有时间维度。当窗口小到单个采样点时,分析也就退化为时域分析,无法做频域分析。

在这里插入图片描述

4 提高频率分辨率的方法

matlab中的spectrogram函数自带两种提高频率分辨率的方法。

第一种是利用功率谱,将能量集中的原理,提高频率分辨率。
用法为:
spectrogram(x,window,noverlap,f,fs,‘reassigned’,‘yaxis’);
前面的同第2节例子,在后面增加了’reassigned’选项。但是这个选项似乎不能用作输出,只能单独以图像形式显示,不知道为什么。

下图为窗函数为256时的两个对比。可以看到右边图的频率分辨率相比左图有了较大提升。但是存在许多噪点。
在这里插入图片描述

第二种是将小于某个数值的数归零,提高频率分辨率。
格式如下:
spectrogram(___,‘MinThreshold’,thresh)
参考第2节中的格式如下
spectrogram(z2,h,wlen-hop,f,fs,‘MinThreshold’,-14,‘yaxis’);
set(gca,‘YDir’,‘reverse’);
结果如下图所示
在这里插入图片描述

5 提取相位信号

短时傅里叶变换不仅可以提取信号的幅值,还可以提取信号的相位。
幅值的提取方法为abs,相位的提取方法为angle。

例:

clear

load mtlb;
z=mtlb;
L=length(z);
t=(1:L)/Fs;

z2=wextend(1,'sym',z,round(length(z)/2));
wlen=512;%设置窗长
hop=1;%平移步长
z2=wkeep1(z2,L+1*wlen);%中间截断


h=hamming(wlen);%设置海明窗的窗长
f=100:1:1000;
[tfr2,f,t2]=spectrogram(z2,h,wlen-hop,f,Fs,'reassigned','yaxis');
tfr2=tfr2*2/wlen*2;
figure
imagesc(t,f,angle(tfr2))%提取相位
colormap([linspace(0,1,32)'*ones(1,3);linspace(1,0,32)'*ones(1,3)])

信号的相位图如下:
在这里插入图片描述

版权声明:本文为博主原创文章,遵循 CC 4.0 BY-SA 版权协议,转载请附上原文出处链接和本声明。
本文链接:https://blog.csdn.net/weixin_42943114/article/details/88735799

智能推荐

使用nginx解决浏览器跨域问题_nginx不停的xhr-程序员宅基地

文章浏览阅读1k次。通过使用ajax方法跨域请求是浏览器所不允许的,浏览器出于安全考虑是禁止的。警告信息如下:不过jQuery对跨域问题也有解决方案,使用jsonp的方式解决,方法如下:$.ajax({ async:false, url: 'http://www.mysite.com/demo.do', // 跨域URL ty..._nginx不停的xhr

在 Oracle 中配置 extproc 以访问 ST_Geometry-程序员宅基地

文章浏览阅读2k次。关于在 Oracle 中配置 extproc 以访问 ST_Geometry,也就是我们所说的 使用空间SQL 的方法,官方文档链接如下。http://desktop.arcgis.com/zh-cn/arcmap/latest/manage-data/gdbs-in-oracle/configure-oracle-extproc.htm其实简单总结一下,主要就分为以下几个步骤。..._extproc

Linux C++ gbk转为utf-8_linux c++ gbk->utf8-程序员宅基地

文章浏览阅读1.5w次。linux下没有上面的两个函数,需要使用函数 mbstowcs和wcstombsmbstowcs将多字节编码转换为宽字节编码wcstombs将宽字节编码转换为多字节编码这两个函数,转换过程中受到系统编码类型的影响,需要通过设置来设定转换前和转换后的编码类型。通过函数setlocale进行系统编码的设置。linux下输入命名locale -a查看系统支持的编码_linux c++ gbk->utf8

IMP-00009: 导出文件异常结束-程序员宅基地

文章浏览阅读750次。今天准备从生产库向测试库进行数据导入,结果在imp导入的时候遇到“ IMP-00009:导出文件异常结束” 错误,google一下,发现可能有如下原因导致imp的数据太大,没有写buffer和commit两个数据库字符集不同从低版本exp的dmp文件,向高版本imp导出的dmp文件出错传输dmp文件时,文件损坏解决办法:imp时指定..._imp-00009导出文件异常结束

python程序员需要深入掌握的技能_Python用数据说明程序员需要掌握的技能-程序员宅基地

文章浏览阅读143次。当下是一个大数据的时代,各个行业都离不开数据的支持。因此,网络爬虫就应运而生。网络爬虫当下最为火热的是Python,Python开发爬虫相对简单,而且功能库相当完善,力压众多开发语言。本次教程我们爬取前程无忧的招聘信息来分析Python程序员需要掌握那些编程技术。首先在谷歌浏览器打开前程无忧的首页,按F12打开浏览器的开发者工具。浏览器开发者工具是用于捕捉网站的请求信息,通过分析请求信息可以了解请..._初级python程序员能力要求

Spring @Service生成bean名称的规则(当类的名字是以两个或以上的大写字母开头的话,bean的名字会与类名保持一致)_@service beanname-程序员宅基地

文章浏览阅读7.6k次,点赞2次,收藏6次。@Service标注的bean,类名:ABDemoService查看源码后发现,原来是经过一个特殊处理:当类的名字是以两个或以上的大写字母开头的话,bean的名字会与类名保持一致public class AnnotationBeanNameGenerator implements BeanNameGenerator { private static final String C..._@service beanname

随便推点

二叉树的各种创建方法_二叉树的建立-程序员宅基地

文章浏览阅读6.9w次,点赞73次,收藏463次。1.前序创建#include&lt;stdio.h&gt;#include&lt;string.h&gt;#include&lt;stdlib.h&gt;#include&lt;malloc.h&gt;#include&lt;iostream&gt;#include&lt;stack&gt;#include&lt;queue&gt;using namespace std;typed_二叉树的建立

解决asp.net导出excel时中文文件名乱码_asp.net utf8 导出中文字符乱码-程序员宅基地

文章浏览阅读7.1k次。在Asp.net上使用Excel导出功能,如果文件名出现中文,便会以乱码视之。 解决方法: fileName = HttpUtility.UrlEncode(fileName, System.Text.Encoding.UTF8);_asp.net utf8 导出中文字符乱码

笔记-编译原理-实验一-词法分析器设计_对pl/0作以下修改扩充。增加单词-程序员宅基地

文章浏览阅读2.1k次,点赞4次,收藏23次。第一次实验 词法分析实验报告设计思想词法分析的主要任务是根据文法的词汇表以及对应约定的编码进行一定的识别,找出文件中所有的合法的单词,并给出一定的信息作为最后的结果,用于后续语法分析程序的使用;本实验针对 PL/0 语言 的文法、词汇表编写一个词法分析程序,对于每个单词根据词汇表输出: (单词种类, 单词的值) 二元对。词汇表:种别编码单词符号助记符0beginb..._对pl/0作以下修改扩充。增加单词

android adb shell 权限,android adb shell权限被拒绝-程序员宅基地

文章浏览阅读773次。我在使用adb.exe时遇到了麻烦.我想使用与bash相同的adb.exe shell提示符,所以我决定更改默认的bash二进制文件(当然二进制文件是交叉编译的,一切都很完美)更改bash二进制文件遵循以下顺序> adb remount> adb push bash / system / bin /> adb shell> cd / system / bin> chm..._adb shell mv 权限

投影仪-相机标定_相机-投影仪标定-程序员宅基地

文章浏览阅读6.8k次,点赞12次,收藏125次。1. 单目相机标定引言相机标定已经研究多年,标定的算法可以分为基于摄影测量的标定和自标定。其中,应用最为广泛的还是张正友标定法。这是一种简单灵活、高鲁棒性、低成本的相机标定算法。仅需要一台相机和一块平面标定板构建相机标定系统,在标定过程中,相机拍摄多个角度下(至少两个角度,推荐10~20个角度)的标定板图像(相机和标定板都可以移动),即可对相机的内外参数进行标定。下面介绍张氏标定法(以下也这么称呼)的原理。原理相机模型和单应矩阵相机标定,就是对相机的内外参数进行计算的过程,从而得到物体到图像的投影_相机-投影仪标定

Wayland架构、渲染、硬件支持-程序员宅基地

文章浏览阅读2.2k次。文章目录Wayland 架构Wayland 渲染Wayland的 硬件支持简 述: 翻译一篇关于和 wayland 有关的技术文章, 其英文标题为Wayland Architecture .Wayland 架构若是想要更好的理解 Wayland 架构及其与 X (X11 or X Window System) 结构;一种很好的方法是将事件从输入设备就开始跟踪, 查看期间所有的屏幕上出现的变化。这就是我们现在对 X 的理解。 内核是从一个输入设备中获取一个事件,并通过 evdev 输入_wayland

推荐文章

热门文章

相关标签