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)])
 
信号的相位图如下：
 

这篇文章首先介绍点的概念，从简单的点到复杂的点，讲解如何用基表示这些点。然后再从复杂的点到函数，引出傅里叶分析，并介绍了傅里叶变换的优缺点，根据其缺点提出改进措施。
 

 
文章目录
 
一、点的概念
1.1 一个简单的点
1.2 基
1.3 函数
   
二、傅里叶分析
三、短时傅里叶变换STFT
四、小波变换
4.1 小波变换的应用
  
 

 
一、点的概念
 
1.1 一个简单的点
 
点在一个直线上可以表示为一个数，在二维平面上可以表示为$x+jy=\alpha e^{j\phi }$。在实数的n维空间中，可以表示为$[x_1,x_2,\dots ,x_n{]}^T$。当有很多点的时候，用上述方式就不太可行了。我们希望用少量的资源，表示出所有的点。那么，如何做呢？
 
1.2 基
 
大家应该都听过基这个名字。那为什么要有基的概念呢？
 
 
 
基：就是用更少的资源，表示全部。
 平凡基（一般坐标系）：$p_1=[1,0,\dots ,0{]}^T,p_2=[0,1,\dots ,0{]}^T,\cdots ,p_n=[0,0,\dots ,1{]}^T$ 空间中任意点，都可以用这个坐标表示。
 
 
举个栗子，现有一点的坐标是$[1,2,3{]}^T$，这个就是又上面$n=3$的平凡基构成的。如果，我想要用更少的资源去表示这个点，我只需要将平凡基换成$p_1=\frac{1}{\sqrt{14}}[1,2,3{]}^T,p_2=\frac{1}{\sqrt{5}}[-2,1,0{]}^T,p_3=\frac{1}{\sqrt{47}}[3,6,-5{]}^T$（其中$p_1,p_2,p_3$是正交的,改变的基不唯一哦）。因此，在这个新基下，同一空间中的点可以表示为$[\sqrt{14},0,0{]}^T$。这样，我们就将3个数表示的点变成了1个数。看下图，可以更加直观的感受。
 
 （蓝色是平凡基，红色是修改后的基）
 
1.3 函数
 
这里给出另一个角度来理解函数。将上面介绍的有限维的点$p(x_1,x_2,\dots ,x_n)$左右延伸，变成$p(\dots ,x_{-1},x_0,x_1,\dots )$无限维，这时无数点联合起来就变成了一个函数。函数又分为周期函数和非周期函数。
 
二、傅里叶分析
 
从函数的角度出发，如何用上面谈到的基来表示任意复杂函数呢？Fourier在1807年提出傅里叶级数。$f(t)={\sum }_{k=-\infty }^{+\infty }{\alpha }_k{e}^{ikt}$，只有确定系数${\alpha }_k$，该函数就可以表达出来了。
 意义：
 1）不同函数的差异就体现在系数${\alpha }_k$上。
 2）对不同特点的函数的分析，可以选择不同的基。
 缺点：
 对于随着时间变换的非平稳信号，它没法区分频谱。只能获取一段信号总体上包含哪些频率的成分，但是对各成分出现的时刻并无所知。因此时域相差很大的两个信号，可能频谱图一样。
 解决办法：
 所以，这里提出两种解决办法。①加入局部参数，比如加窗。②换基。
 关于傅里叶级数和傅里叶变换的具体公式推导总结，见https://blog.csdn.net/weixin_46017950/article/details/114691667。
 
三、短时傅里叶变换STFT
 
方法一：加窗。
 STFT在傅里叶变换的基础上加窗，分段做FFT变换，假定认为信号在窗宽度的时间内是平稳的。但是，窗宽度太窄，信号太短信息少，频率分析不准确，频率分辨率差；窗太宽时间分辨率差。对于时变的非稳态信号，高频适合小窗口，低频适合大窗口。（高频，顾名思义就是信号变化很快，所以，在时域中，需要小窗口，时间分辨率高。）然而STFT的窗口是固定的，在一次STFT中宽度不会变化，所以STFT还是无法满足非稳态信号变化的频率的需求。不采用可变窗的STFT，是因为这样做冗余会太严重，STFT做不到正交化，这也是它的一大缺陷。
 
四、小波变换
 
方法二：换基。
 何为小波呢？“小”是指衰减性，比如有些小波基只有局部是非零，这也称为紧支性。“波”是指波动性。
 小波变换直接把傅里叶变换的基给换了——将无限长的三角函数基换成了有限长的会衰减的小波基。这样不仅能够获取频率，还可以定位到时间，获取时频分析。目标：时间分辨率和频率分辨率可以随着本身信号特变进行自适应调整。这里的“自适应”是很难的。
 $$WT(a,b)=\frac{1}{\sqrt{a}}{\int }_{-\infty }^{\infty }f(t)\ast \psi \left(\frac{t-b}{a}\right)dt$$
 小波还有一些好处，比如，我们知道对于突变信号，傅里叶变换存在吉布斯效应，我们用有限长的频率信号是怎么也拟合不好突变信号的。
 
JPEG2000压缩就是用正交小波变换。比如典型的正交基：二维笛卡尔坐标系的（1,0）、（0,1），用它们表达一个信号显然非常高效，计算简单。而如果用三个互成120°的向量表达，则会有信息冗余，有重复表达。但是并不意味着正交一定优于不正交。比如，如果是做图像增强，有时候反而希望能有一些冗余信息，更利于对噪声的抑制和对某些特征的增强。
 
4.1 小波变换的应用
 
1、数据压缩。目前许多应用领域（如卫星监测、地震勘探、天气预报）都存在海量数据传输或存储问题，如果不对数据进行压缩，数量巨大的数据就很难存储、处理和传输。因此，伴随小波分析的诞生，数据压缩一直是小波分析的重要应用领域之一，并由此带来巨大的经济效益和社会效益。
 2、语音分析与处理。小波理论应用于语音分析与处理的主要内容包括：清/浊音分割；基音检测与声门开启时刻定位；去噪、压缩、重建几个方面。
 3、瞬态信号或图像的突变点常包含有很重要的故障信息，例如，机械故障、电力系统故障、脑电图、心电图中的异常、地下目标的位置及形状等，都对应于测试信号的突变点。因此，小波分析在故障检测和信号的多尺度边缘特征提取方面的应用具有广泛的应用前景。
 4、神经网络与小波分析相结合，分形几何与小波分析相结合是国际上研究的热点之一。基于神经网络的智能处理技术，模糊计算、进化计算与神经网络结合的研究，没有小波理论的嵌入很难取得突破。非线性科学的研究正呼唤小波分析，也许非线性小波分析是解决非线性科学问题的理性工具。
 
下一篇文章写一写小波变换在图像方面的分析与应用（Matlab代码）。

 
 


 
作者丨咚懂咚懂咚@知乎（已授权）
 
来源丨https://zhuanlan.zhihu.com/p/22450818
 
转载丨极市平台
 
导读
 
 
想要正确的认识小波变换就必须先了解傅里叶变换，本文作者按照傅里叶-短时傅里叶变换-小波变换的顺序，由浅到深的解释小波变换的缘由以及思路。帮助初学者们深入理解傅里叶变换和小波变换。
 
从傅里叶变换到小波变换，并不是一个完全抽象的东西，可以讲得很形象。小波变换有着明确的物理意义，如果我们从它的提出时所面对的问题看起，可以整理出非常清晰的思路。
 
下面我就按照傅里叶-->短时傅里叶变换-->小波变换的顺序，讲一下为什么会出现小波这个东西、小波究竟是怎样的思路。（反正题主要求的是通俗形象，没说简短，希望不会太长不看。。）
 
一、傅里叶变换
 
关于傅里叶变换的基本概念在此我就不再赘述了，默认大家现在正处在理解了傅里叶但还没理解小波的道路上。（在第三节小波变换的地方我会再形象地讲一下傅里叶变换）
 
下面我们主要讲傅里叶变换的不足。即我们知道傅里叶变化可以分析信号的频谱，那么为什么还要提出小波变换？答案就是@方沁园所说的，“对非平稳过程，傅里叶变换有局限性”。看如下一个简单的信号：
 
 
做完FFT（快速傅里叶变换）后，可以在频谱上看到清晰的四条线，信号包含四个频率成分。
 
一切没有问题。但是，如果是频率随着时间变化的非平稳信号呢？
 
 
如上图，最上边的是频率始终不变的平稳信号。而下边两个则是频率随着时间改变的非平稳信号，它们同样包含和最上信号相同频率的四个成分。
 
做FFT后，我们发现这三个时域上有巨大差异的信号，频谱（幅值谱）却非常一致。尤其是下边两个非平稳信号，我们从频谱上无法区分它们，因为它们包含的四个频率的信号的成分确实是一样的，只是出现的先后顺序不同。
 
可见，傅里叶变换处理非平稳信号有天生缺陷。它只能获取一段信号总体上包含哪些频率的成分，但是对各成分出现的时刻并无所知。因此时域相差很大的两个信号，可能频谱图一样。
 
然而平稳信号大多是人为制造出来的，自然界的大量信号几乎都是非平稳的，所以在比如生物医学信号分析等领域的论文中，基本看不到单纯傅里叶变换这样naive的方法。
 
 
上图所示的是一个正常人的事件相关电位。对于这样的非平稳信号，只知道包含哪些频率成分是不够的，我们还想知道各个成分出现的时间。知道信号频率随时间变化的情况，各个时刻的瞬时频率及其幅值——这也就是时频分析。
 
二、短时傅里叶变换（Short-time Fourier Transform, STFT）
 
一个简单可行的方法就是——加窗。我又要套用方沁园同学的描述了，“把整个时域过程分解成无数个等长的小过程，每个小过程近似平稳，再傅里叶变换，就知道在哪个时间点上出现了什么频率了。”这就是短时傅里叶变换。
 
看图：
 
 
时域上分成一段一段做FFT，不就知道频率成分随着时间的变化情况了吗！
 
用这样的方法，可以得到一个信号的时频图了：
 
 
——此图像来源于“THE WAVELET TUTORIAL”
 
图上既能看到10Hz, 25 Hz, 50 Hz, 100 Hz四个频域成分，还能看到出现的时间。两排峰是对称的，所以大家只用看一排就行了。
 
是不是棒棒的？时频分析结果到手。但是STFT依然有缺陷。
 
使用STFT存在一个问题，我们应该用多宽的窗函数？
 
窗太宽太窄都有问题：
 
 
 
窗太窄，窗内的信号太短，会导致频率分析不够精准，频率分辨率差。窗太宽，时域上又不够精细，时间分辨率低。
 
（这里插一句，这个道理可以用海森堡不确定性原理来解释。类似于我们不能同时获取一个粒子的动量和位置，我们也不能同时获取信号绝对精准的时刻和频率。这也是一对不可兼得的矛盾体。我们不知道在某个瞬间哪个频率分量存在，我们知道的只能是在一个时间段内某个频带的分量存在。所以绝对意义的瞬时频率是不存在的。）
 
看看实例效果吧：
 
 
 
 
——此图像来源于“THE WAVELET TUTORIAL”
 
上图对同一个信号（4个频率成分）采用不同宽度的窗做STFT，结果如右图。用窄窗，时频图在时间轴上分辨率很高，几个峰基本成矩形，而用宽窗则变成了绵延的矮山。但是频率轴上，窄窗明显不如下边两个宽窗精确。
 
所以窄窗口时间分辨率高、频率分辨率低，宽窗口时间分辨率低、频率分辨率高。对于时变的非稳态信号，高频适合小窗口，低频适合大窗口。然而STFT的窗口是固定的，在一次STFT中宽度不会变化，所以STFT还是无法满足非稳态信号变化的频率的需求。
 
三、小波变换
 
那么你可能会想到，让窗口大小变起来，多做几次STFT不就可以了吗？！没错，小波变换就有着这样的思路。
 
但事实上小波并不是这么做的（关于这一点，方沁园同学的表述“小波变换就是根据算法，加不等长的窗，对每一小部分进行傅里叶变换”就不准确了。小波变换并没有采用窗的思想，更没有做傅里叶变换。）
至于为什么不采用可变窗的STFT呢，我认为是因为这样做冗余会太严重，STFT做不到正交化，这也是它的一大缺陷。
 
于是小波变换的出发点和STFT还是不同的。STFT是给信号加窗，分段做FFT；而小波直接把傅里叶变换的基给换了——将无限长的三角函数基换成了有限长的会衰减的小波基。这样不仅能够获取频率，还可以定位到时间了~
 
【解释】
 
来我们再回顾一下傅里叶变换吧，没弄清傅里叶变换为什么能得到信号各个频率成分的同学也可以再借我的图理解一下。
 
傅里叶变换把无限长的三角函数作为基函数：
 
 
这个基函数会伸缩、会平移（其实本质并非平移，而是两个正交基的分解）。缩得窄，对应高频；伸得宽，对应低频。然后这个基函数不断和信号做相乘。某一个尺度（宽窄）下乘出来的结果，就可以理解成信号所包含的当前尺度对应频率成分有多少。于是，基函数会在某些尺度下，与信号相乘得到一个很大的值，因为此时二者有一种重合关系。那么我们就知道信号包含该频率的成分的多少。
 
仔细体会可以发现，这一步其实是在计算信号和三角函数的相关性。
 
 
 
看，这两种尺度能乘出一个大的值（相关度高），所以信号包含较多的这两个频率成分，在频谱上这两个频率会出现两个峰。
 
以上，就是粗浅意义上傅里叶变换的原理。
 
如前边所说，小波做的改变就在于，将无限长的三角函数基换成了有限长的会衰减的小波基。
 
 
这就是为什么它叫“小波”，因为是很小的一个波嘛~
 
 
从公式可以看出，不同于傅里叶变换，变量只有频率ω，小波变换有两个变量：尺度a（scale）和平移量 τ（translation）。尺度a控制小波函数的伸缩，平移量 τ控制小波函数的平移。尺度就对应于频率（反比），平移量 τ就对应于时间。
 
 
当伸缩、平移到这么一种重合情况时，也会相乘得到一个大的值。这时候和傅里叶变换不同的是，这不仅可以知道信号有这样频率的成分，而且知道它在时域上存在的具体位置。
 
而当我们在每个尺度下都平移着和信号乘过一遍后，我们就知道信号在每个位置都包含哪些频率成分。
 
看到了吗？有了小波，我们从此再也不害怕非稳定信号啦！从此可以做时频分析啦！
 
做傅里叶变换只能得到一个频谱，做小波变换却可以得到一个时频谱！
 
 
↑：时域信号
 
 
↑：傅里叶变换结果
 
 
——此图像来源于“THE WAVELET TUTORIAL”
↑：小波变换结果
 
小波还有一些好处，比如，我们知道对于突变信号，傅里叶变换存在吉布斯效应，我们用无限长的三角函数怎么也拟合不好突变信号：
 
 
然而衰减的小波就不一样了：
 
 
以上，就是小波的意义。
 
-----------------------------------------------------------------------------------------------------------
 
以上只是用形象地给大家展示了一下小波的思想，希望能对大家的入门带来一些帮助。毕竟如果对小波一无所知，直接去看那些堆砌公式、照搬论文语言的教材，一定会痛苦不堪。
在这里推荐几篇入门读物，都是以感性介绍为主，易懂但并不深入，对大家初步理解小波会很有帮助。文中有的思路和图也选自于其中：
1. THE WAVELET TUTORIAL （强烈推荐，点击链接：Ihttp://users.rowan.edu/~polikar/WTtutorial.html）
2. WAVELETS：SEEING THE FOREST AND THE TREES
3. A Really Friendly Guide to Wavelets
4. Conceptual wavelets
 
但是真正理解透小波变换，这些还差得很远。比如你至少还要知道有一个“尺度函数”的存在，它是构造“小波函数”的关键，并且是它和小波函数一起才构成了小波多分辨率分析，理解了它才有可能利用小波做一些数字信号处理；你还要理解离散小波变换、正交小波变换、二维小波变换、小波包……这些内容国内教材上讲得也很糟糕，大家就一点一点啃吧~
 
------
 
一些问题的回答：
 
1. 关于海森堡不确定性原理
 
不确定性原理，或者叫测不准原理，最早出自量子力学，意为在微观世界，粒子的位置与动量不可同时被确定。但是这个原理并不局限于量子力学，有很多物理量都有这样的特征，比如能量和时间、角动量和角度。体现在信号领域就是时域和频域。不过更准确一点的表述应该是：一个信号不能在时空域和频域上同时过于集中；一个函数时域越“窄”，它经傅里叶变换的频域后就越“宽”。
 
如果有兴趣深入研究一下的话，这个原理其实非常耐人寻味。信号处理中的一些新理论在根本上也和它有所相连，比如压缩感知。如果你剥开它复杂的数学描述，最后会发现它在本质上能实现其实和不确定性原理密切相关。而且大家不觉得这样一些矛盾的东西在哲学意义上也很奇妙吗？
 
2. 关于正交化
 
什么是正交化？为什么说小波能实现正交化是优势?
 
简单说，如果采用正交基，变换域系数会没有冗余信息，变换前后的信号能量相等，等于是用最少的数据表达最大的信息量，利于数值压缩等领域。JPEG2000压缩就是用正交小波变换。
 
比如典型的正交基：二维笛卡尔坐标系的（1,0）、（0,1），用它们表达一个信号显然非常高效，计算简单。而如果用三个互成120°的向量表达，则会有信息冗余，有重复表达。
 
但是并不意味着正交一定优于不正交。比如如果是做图像增强，有时候反而希望能有一些冗余信息，更利于对噪声的抑制和对某些特征的增强。
 
3. 关于瞬时频率
 
原问题：图中时刻点对应一频率值，一个时刻点只有一个信号值，又怎么能得到他的频率呢？
 
很好的问题。如文中所说，绝对意义的瞬时频率其实是不存在的。单看一个时刻点的一个信号值，当然得不到它的频率。我们只不过是用很短的一段信号的频率作为该时刻的频率，所以我们得到的只是时间分辨率有限的近似分析结果。这一想法在STFT上体现得很明显。小波用衰减的基函数去测定信号的瞬时频率，思想也类似。（不过到了Hilbert变换，思路就不一样了，以后有机会细讲）
 
4. 关于小波变换的不足
 
这要看和谁比了。
 
A.作为图像处理方法，和多尺度几何分析方法（超小波）比：
对于图像这种二维信号的话，二维小波变换只能沿2个方向进行，对图像中点的信息表达还可以，但是对线就比较差。而图像中最重要的信息恰是那些边缘线，这时候ridgelet（脊波）, curvelet（曲波）等多尺度几何分析方法就更有优势了。
 
B. 作为时频分析方法，和希尔伯特-黄变换（HHT）比：
相比于HHT等时频分析方法，小波依然没脱离海森堡测不准原理的束缚，某种尺度下，不能在时间和频率上同时具有很高的精度；以及小波是非适应性的，基函数选定了就不改了。
 
5. 关于文中表述的严谨性
 
评论中有不少朋友提到，我的一些表述不够精准。这是肯定的，并且我也是知道的。比如傅里叶变换的理解部分，我所说的那种“乘出一个大的值”的表述肯定是不够严谨的。具体我也在评论的回答中做了解释。我想说的是通俗易懂和精确严谨实在难以兼得，如果要追求严谨，最好的就是教科书上的数学表达，它们无懈可击，但是对于初学者来说，恐怕存在门槛。如果要通俗解释，必然只能侧重一个关键点，而出现漏洞。我想这也是教科书从来不把这些通俗解释写出来的原因吧——作者们不是不懂，而是怕写错。所以想深入理解傅里叶变换和小波变换的朋友还请认真学习教材，如果这篇文章能给一些初学者一点点帮助，我就心满意足了。