短时傅里叶变换（5）

目录：

前言

之前讲了时频分析的原理，现在来讲讲它在matlab里面的实现。
想要复习原理的同学，可以参照一下这篇：短时傅里叶分析
本次讲解终阶的函数使用，
基础的可以参见前面的：短时傅里叶实现（1）
中阶的可以参见前面的：短时傅里叶实现（2）
高阶的可以参见前面的：短时傅里叶实现（3）
进阶的可以参见前面的：短时傅里叶实现（4）

实验环境

本文的所有实验都是在matlab2016a下通过的。

Matlab spectrogram函数

谱图函数：使用短时傅里叶变换化成短时傅里叶变换的谱图。

1语法

标准函数引用方式如下：

[___] = spectrogram(___,freqrange) 
[___] = spectrogram(___,spectrumtype) 
[___] = spectrogram(___,'MinThreshold',thresh) 
 
 

 
 
1
2
3

2使用说明

[] = spectrogram(,freqrange)
返回由frequencge确定范围的功率谱密度或者功率谱估计，正确的选项有单边（oneside），双边（twoside），中心（centered）。
[] = spectrogram(,spectrumtype)
返回功率谱估计。如果使用psd参数。返回能量，如果使用power参数
[] = spectrogram(,’MinThreshold’,thresh)
将门限值之以下的信号，ps 10log10log10 thresh .设置为0.

3代码如下

3.1谱图聚集和门限设置

%产生一个信号，采样频率是1khz，采样时间是两秒                                                                                                                 
nSamp=2048;
Fs=1024;
t=(0:nSamp-1)'/Fs;
%在第一秒，信号包含一个400hz的正弦和凹二次鸟声信号，鸟声信号关于间隔中点是对称的，起始和结束为250hz，
%包%含一个最小的150hz信号
t1=t(1:nSamp/2);
x11=sin(2*pi*400*t1);
x12=chirp(t1-t1(nSamp/4),150,nSamp/Fs,1750,'quadratic');
x1=x11+x12;
% 信号的剩余部分包含两个线性衰落鸟声成分，一个鸟声的初始频率是250hz最后降到100hz，其他线性鸟声有一个出声频率400hz，并减到250hz。
t2=t(nSamp/2+1:nSamp);
x21=chirp(t2,400,nSamp/Fs,100);
x22=chirp(t2,550,nSamp/Fs,250);
x2=x21+x22;
%给信号添加高斯白噪声，信噪比为20db，设置随机数生成器，为了可重复的结果。
SNR=20；
rng('default')
sig=[x1;x2];
sig=sig+randn(size(sig))*std(sig)/db2mag(SNR);
%计算并绘制信号的谱图，长度为63的kaizer窗，形状参数17，10个左右重叠值，FFT长度为256.
nwin=63;
wind=kaiser(nwin,17);
nlap=nwin-10;
nfft=256;
spectrogram(sig,wind,nlap,nfft,Fs,'yaxis')
 
 

 
 
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26

如图所示

%任何能量值小于噪声（SNR）的信号都被设置为0

spectrogram(sig,wind,nlap,nfft,Fs,’MinThreshold’,-SNR,’yaxis’)

%将谱功率密度估计值压缩到能量中心的位置
spectrogram(sig,wind,nlap,nfft,Fs,’reassign’,’yaxis’)
结果如图所示：

%对集中能量谱的信号设置阈值，每个元素能量值低于SNR的被设置为0
spectrogram(sig,wind,nlap,nfft,Fs,’reassign’,’MinThreshold’,-SNR,’yaxis’)
如图所示：

函数相关功能展示完毕，大家可以根据自己需要，选择参数

---

参考：

mathworks