对音频信号作短时傅里叶变换（STFT）处理，并绘制语谱图
摘要：录制一段音频，分别在matlab，python两种环境下，对其作短时傅里叶变换（STFT），最终得到期望的语谱图。
一、前言
1. 基础概念：
在对音频信号进行分析处理前，先简要回顾一下所用到的分析函数傅里叶变换的相关知识。
什么是傅里叶变换？
傅里叶的基本定义与性质在这里就不作赘述了，文章主要想说明它的主要应用，以助于大家对这个概念有一个更为形象的认识。（傅里叶基础概念与性质，推荐观看b站李永乐的讲解，文末附视频链接）
傅里叶变换是一种分析信号的方法，它可分析信号的成分，也可用这些成分合成信号。在分析信号时，主要应用于处理平稳信号，通过傅里叶变换可以获取一段信号总体上包含哪些频率的成分，但是对各成分出现的时刻无法得知。
因此对于非平稳信号，傅里叶变换就显示出了它的局限性，而我们日常生活中的绝大多数音频都是非平稳信号的。而解决这一问题的方法，就是采用短时傅里叶变换或者小波变换，对信号进行处理。
什么是短时傅里叶变换？
短时傅里叶变换（STFT）的核心思想：“加窗”，即把整个时域过程分解成无数个等长的小过程，每个小过程近似平稳，再对每个小过程进行傅里叶变换（FFT）。尽管STFT可以处理非平稳信号，但是它仍然有其局限性，即对窗函数的宽窄无法做到精确定义。
窗函数选择太窄，窗内的信号太短，会导致频率的分析不够精准，频率分辨率差；窗选的太宽，时域上又不够精细，时间分辨率低
傅里叶与短时傅里叶的联系与区别？
什么是小波变换？
小波变换的核心思想：“把傅里叶变换的无限长三角函数的基换成有限长的会衰减的小波基”，这样不仅可以获取频率，还可以定位时间。更为详细的介绍可以参考本段末附上的参考来源，查询大佬相关阐述。
本文参考来源：
/article/id-xiaobobianhuan?p=1 傅里叶–短时傅里叶–小波分析
/daaikuaichuan/article/details/80781505 小波与短时傅里叶的区别
/view/4b9bb22c30b765ce0508763231126edb6f1a768e.html?fr=search-1_income7 短时傅里叶的概念理解
什么是语谱图？
语谱图：时间依赖于傅里叶分析的显示图形，实际上是一种动态频谱，综合了频谱图和时域波形图的优点，明显地显示出语音频谱随时间的变化情况。（其中，纵轴为频率，横着为时间，任一给定频率成分在给定时刻的强弱，用点的黑白度来表示）
2. 概念的形象化理解：
短时傅里叶变换进行音频分析流程的直观理解
欧拉公式可以将任意函数转化为一系列正（余）弦函数之和。任意函数在这里指我们的时域信号，而正（余）弦函数包含信号的频率和对应振幅信息。（时域处理）傅立叶变换可以将时间0～t内采集的信号（时域，横轴时间、纵轴大小）分解为不同频率上的信号分量（频域，横轴频率、纵轴大小）。短时傅立叶采用滑动窗口机制，设定窗口大小和步长，让窗口在时域信号上滑动，分别计算每个窗口的傅立叶变换，形成了不同时间窗口对应的频域信号，拼接起来就成为了频率随时间变化的数据（时频信号，数据绘图得语谱图）
文章参考来源：
/p/7e160442830f 直观理解短时傅里叶变换 STFT（简书）
/video/BV1A4411Y7vj/?spm_id_from=333.788.videocard.1 傅里叶变换-李永乐（b站视频，讲的很好）
3. 思路流程：
编程思路：
短时傅里叶变换（STFT）就是先把一个函数和窗函数进行相乘，然后再进行一维的傅里叶变换。并通过窗函数的滑动得到一系列的傅里叶变换结果，将这些结果排开便得到一个二维的表象。具体而言，可分为如下几步（以matlab源码为例）：
（1）读取音频文件。（调用audioread，处理得到一个保存音频数据的数组，和一个采样频率）
（2）确定相关参数。（如窗函数、窗长、重叠点数，重叠长度，傅里叶点数等）
（3）调用spectrogram函数，做短时傅里叶变换。（S-将输入信号做STFT处理后得到的二维含时间、频率序列的数组数据；更为详细的处理过程见spectrogram函数定义）
（4）根据处理后的时频矩阵，绘制语谱图.
文章参考来源：
/zhaoyinhui0802/article/details/53048362 短时傅里叶变换原理解
整体思路：
录制一段音频或寻找一段音频（.MP3/.m4a）,通过文件格式转换为.wav文件（网上有在线转换的方式，文后会给出网址链接），在python（调用librosa函数）或者matlab（调用audioread，spectogram函数）中对文件进行处理（STFT变换），最后对处理后的信号作绘图输出，得到语谱图。
文章参考来源：
https://cn.office-/ 文件在线转换网址
4. 必要准备：
音频文件+python应用环境+matlab应用环境
python应用环境：
作者是在vscode中编写python代码，CSDN中有许多博主介绍了vscode使用python的方法，在这里不做多的赘述，文章用到的调用函数为librosa，文章后面会贴出librosa库安装的说明。 可通过pip list 或 conda list 查询python 或anaconda是否包含librosa 库文件
文章参考来源：
/Small_Yogurt/article/details/104964760 vscode安装教程
/qq_40197828/article/details/93468787 vscode中配置python环境（可以选择配置anaconda中的python，这样会减少后续代码实现阶段，函数库的安装，如numpy,librosa等）
/u010824101/article/details/85239223 Windows10中用Anaconda的conda环境和VSCode工具来编写python代码
/zzc15806/article/details/79603994 音频处理库—librosa的安装与使用
matlab应用环境：
如何安装以及破解matlab，请参考下方资料来源（大佬写的非常详细），这里不作赘述了。
安装前须知（提醒）：
1.安装全程须断开网络，否则破解不成功；
2.解压和安装前先关闭360、腾讯管家等杀毒软件，防止误杀破解补丁，导致破解失败；
3.Matlab b适用于WIN7/8/10（64位）系统，亲测可用！
4.推荐电脑最低配置：内存8G+，处理器酷睿I5+；
文章参考来源：
/qq_26900233/article/details/88816789 MATLAB Rb 安装教程（含软件资源）
二、代码实现
python源码：import librosaimport librosa.displayimport matplotlib.pyplot as plt# 读取音频文件filepath = 'C:\\Users\\非黑不即白\\Desktop\\Audio_processing_STFT\\python_project\\data\\'filename = filepath + 'test_digital.wav'x, sr = librosa.load(filename, sr=None) # x--音频时间序列(一维数组) ； sr--音频的采样率# STFT处理绘制声谱图X = librosa.stft(x)Xdb = librosa.amplitude_to_db(abs(X)) # X--二维数组数据plt.figure(figsize=(14, 5))librosa.display.specshow(Xdb, sr=sr, x_axis='time', y_axis='log')plt.colorbar()plt.title('STFT transform processing audio signal')plt.show()
matlab源码：
[X, Fs]=audioread('C:\Users\非黑不即白\Desktop\Audio_processing_STFT\matlab_project\test_digital.wav'); % Fs 采样率 48000% audioread函数读取音频文件（X--保存音频信号的数据；Fs--音频采样率）figure; % 绘图win_sz = 128; % 窗函数长度设置为128han_win = hanning(win_sz); % 选择海明窗nfft = 2^nextpow2(length(han_win)); % DFT点数（通常取最接近信号长度的2的整数次幂，即 nfft = 2^nextpow2(length(window)) ）=128=win_sizenooverlap = win_sz - 1; % 重叠长度（而随着noverlap参数的引入，增大了时间轴分辨率，即每隔（Nw - noverlap）长度进行一次频率轴的更新，随着noverlap逐渐接近Nw，图像上表现为时间轴更加“细腻”，但随之而来的肯定是计算次数的增加。）[S, F, T, P] = spectrogram(X(:,1), window, nooverlap, nfft, Fs); % spectrogram注释：x：输入的信号 向量；window：窗口长度 该函数默认使用海明窗；noverlap：各段之间重叠的采样点数 即两窗口相重叠的部分；nfft：对窗口下的信号做FFT的点数；fs：信号的采样率% s：输入信号x的短时傅里叶变换。（它的每一列包含一个短期局部时间的频率成分估计，时间沿列增加，频率沿行增加。）% F：为四舍五入的频率，其长度等于S的行数。% T：频谱图计算的时刻点，值为窗的时刻中值。% P：功率谱密度PSD(Power Spectral Density) imagesc(T, F, log10(abs(S))) % 绘图函数（将矩阵S中的元素数值按大小转化为不同颜色，并在坐标轴对应位置处以这种颜色染色；T、F为其横纵坐标）colorbar；set(gca, 'YDir', 'normal')% 坐标轴设置（ydir表示Y轴；normal-坐标轴正序；reverse--倒序）xlabel('Time (secs)')ylabel('Freq (Hz)')title('STFT transform spectrum')
小波变换处理（matlab）:
[X, Fs]=audioread('C:\Users\非黑不即白\Desktop\Audio_processing_STFT\matlab_project\test_digital.wav'); % Fs 采样率 48000% audioread函数读取音频文件（X--保存音频信号的数据；Fs--音频采样率）wavename='cmor3-3'; % 选定小波基totalscal=256;Fc=centfrq(wavename); % 小波的中心频率 测得Fc = 3c=2*Fc*totalscal; % 测得c = 1536scals=c./(1:totalscal); % 求得尺度f=scal2frq(scals,wavename,1/Fs); % 将尺度转换为频率 % 频率在0-500Hz取1024<span style="font-family:Arial, Helvetica, sans-serif;">个点</span>coefs = cwt(X(:,1),scals,wavename); % 求连续小波系数t=0:1/Fs:size(X(:,1))/Fs;figureimagesc(t,f,abs(coefs));set(gca,'YDir','normal')colorbar;xlabel('时间 t/s');ylabel('频率 f/Hz');title('小波时频图')
文章参考来源：
/lanchunhui/article/details/72240693 matlab 时频分析（短时傅里叶变换、STFT）
/weixin_42765703/article/details/104871604 解析MATLAB短时傅里叶变换函数spectrogram()
/zb1165048017/article/details/80682473 【音频处理】短时傅里叶变换
/weixin_44586750/article/details/103219952 数字信号处理——时频分析（短时傅里叶变换+手动实现时频分析）
/qq_42815385/article/details/89095135 STFT（短时傅里叶变换）音频特征提取，用于语音识别 python
/Barry_J/article/details/81065932 音频特征提取——python/ librosa工具包使用
/qq_34218078/article/details/101255636 librosa语音信号处理
/yunman/article/details/103542690 librosa处理音频信号
三、总结（作者小白话）
写这篇博客的缘由，主要想对近期音频信号处理相关知识的学习，做一个系统性的梳理，同时也希望能够对大家有所帮助（学习阶段，查阅了不少大佬文章，收益匪浅，这些都成为我写博客，分享自己的学习历程，与经验教训的一个动力）。整个学习阶段，从了解傅里叶变换，到python、matlab环境的配置和使用，以及后续的音频文件处理，STFT变换等，走过了许多弯路，也遇到了一些不常见的bug（主要是python环境配置，库安装时），但最终还是坚持了下来。目前对音频处理的相关学习还不是很深入，文章还存在许多不全面的地方，在后续深入学习过程中，会对不合理的地方做出调整，也希望各位大佬能多多指教。相互学习，共同进步啊，哈哈哈。（文章引用了一些大佬的文章（文章小节处有标注），以便各位同仁能更全面的了解，当然也会存在一些疏漏，请指教，望见谅。）补充，写文章阶段时，发现了一位大佬关于python语音基础操作的文章，内容很详实，包含音频生成，语音分析、处理的介绍，收益匪浅，文章最后附原文链接，希望有助于大家理解语音信号处理。文章参考来源：
/sinat_18131557/article/details/106440757 Python语音基础操作（对python处理语音信号的过程介绍的很详细，推荐！）
