對音頻信號做短時傅里葉變換（STFT）/小波變換處理（python + matlab）

對音頻信號做短時傅里葉變換（STFT）處理，並繪製語譜圖

摘要：錄製一段音頻，分別在matlab，python兩種環境下，對其做短時傅里葉變換（STFT），最終獲得指望的語譜圖。html

1、前言

1. 基礎概念：
在對音頻信號進行分析處理前，先簡要回顧一下所用到的分析函數傅里葉變換的相關知識。
python

什麼是傅里葉變換？
傅里葉的基本定義與性質在這裏就不做贅述了，文章主要想說明它的主要應用，以助於你們對這個概念有一個更爲形象的認識。（傅里葉基礎概念與性質，推薦觀看b站李永樂的講解，文末附視頻連接）

git
傅里葉變換是一種分析信號的方法，它可分析信號的成分，也可用這些成分合成信號。在分析信號時，主要應用於處理平穩信號，經過傅里葉變換能夠獲取一段信號整體上包含哪些頻率的成分，可是對各成分出現的時刻沒法得知。編程
所以對於非平穩信號，傅里葉變換就顯示出了它的侷限性，而咱們平常生活中的絕大多數音頻都是非平穩信號的。而解決這一問題的方法，就是採用短時傅里葉變換或者小波變換，對信號進行處理。數組
什麼是短時傅里葉變換？

短時傅里葉變換（STFT）的核心思想：「加窗」，即把整個時域過程分解成無數個等長的小過程，每一個小過程近似平穩，再對每一個小過程進行傅里葉變換（FFT）。儘管STFT能夠處理非平穩信號，可是它仍然有其侷限性，即對窗函數的寬窄沒法作到精肯定義。

窗函數選擇太窄，窗內的信號過短，會致使頻率的分析不夠精準，頻率分辨率差；窗選的太寬，時域上又不夠精細，時間分辨率低

網絡
傅里葉與短時傅里葉的聯繫與區別？

ide
什麼是小波變換？函數

小波變換的核心思想：「把傅里葉變換的無限長三角函數的基換成有限長的會衰減的小波基」，這樣不只能夠獲取頻率，還能夠定位時間。更爲詳細的介紹能夠參考本段末附上的參考來源，查詢大佬相關闡述。工具

本文參考來源：
https://rf.eefocus.com/article/id-xiaobobianhuan?p=1 傅里葉–短時傅里葉–小波分析
https://wenku.baidu.com/view/4b9bb22c30b765ce0508763231126edb6f1a768e.html?fr=search-1_income7 短時傅里葉的概念理解

學習

什麼是語譜圖？

語譜圖：時間依賴於傅里葉分析的顯示圖形，其實是一種動態頻譜，綜合了頻譜圖和時域波形圖的優勢，明顯地顯示出語音頻譜隨時間的變化狀況。（其中，縱軸爲頻率，橫着爲時間，任一給定頻率成分在給定時刻的強弱，用點的黑白度來表示）

2. 概念的形象化理解：

短時傅里葉變換進行音頻分析流程的直觀理解
歐拉公式能夠將任意函數轉化爲一系列正（餘）弦函數之和。任意函數在這裏指咱們的時域信號，而正（餘）弦函數包含信號的頻率和對應振幅信息。（時域處理）
傅立葉變換能夠將時間0～t內採集的信號（時域，橫軸時間、縱軸大小）分解爲不一樣頻率上的信號份量（頻域，橫軸頻率、縱軸大小）。
短時傅立葉採用滑動窗口機制，設定窗口大小和步長，讓窗口在時域信號上滑動，分別計算每一個窗口的傅立葉變換，造成了不一樣時間窗口對應的頻域信號，拼接起來就成爲了頻率隨時間變化的數據（時頻信號，數據繪圖得語譜圖）

文章參考來源：
https://www.jianshu.com/p/7e160442830f 直觀理解短時傅里葉變換 STFT（簡書）
https://www.bilibili.com/video/BV1A4411Y7vj/?spm_id_from=333.788.videocard.1 傅里葉變換-李永樂（b站視頻，講的很好）

3. 思路流程：

編程思路：

短時傅里葉變換（STFT）就是先把一個函數和窗函數進行相乘，而後再進行一維的傅里葉變換。並經過窗函數的滑動獲得一系列的傅里葉變換結果，將這些結果排開便獲得一個二維的表象。具體而言，可分爲以下幾步（以matlab源碼爲例）：
（1）讀取音頻文件。（調用audioread，處理獲得一個保存音頻數據的數組，和一個採樣頻率）
（2）肯定相關參數。（如窗函數、窗長、重疊點數，重疊長度，傅里葉點數等）
（3）調用spectrogram函數，作短時傅里葉變換。（S-將輸入信號作STFT處理後獲得的二維含時間、頻率序列的數組數據；更爲詳細的處理過程見spectrogram函數定義）
（4）根據處理後的時頻矩陣，繪製語譜圖.

文章參考來源：
https://blog.csdn.net/zhaoyinhui0802/article/details/53048362 短時傅里葉變換原理解

總體思路：

錄製一段音頻或尋找一段音頻（.MP3/.m4a）,經過文件格式轉換爲.wav文件（網上有在線轉換的方式，文後會給出網址連接），在python（調用librosa函數）或者matlab（調用audioread，spectogram函數）中對文件進行處理（STFT變換），最後對處理後的信號做繪圖輸出，獲得語譜圖。

文章參考來源：
https://cn.office-converter.com/ 文件在線轉換網址

4. 必要準備：

音頻文件+python應用環境+matlab應用環境

python應用環境：

做者是在vscode中編寫python代碼，CSDN中有許多博主介紹了vscode使用python的方法，在這裏不作多的贅述，文章用到的調用函數爲librosa，文章後面會貼出librosa庫安裝的說明。

可經過pip list 或 conda list 查詢python 或anaconda是否包含librosa 庫文件

文章參考來源：
https://blog.csdn.net/Small_Yogurt/article/details/104964760 vscode安裝教程
https://blog.csdn.net/qq_40197828/article/details/93468787 vscode中配置python環境（能夠選擇配置anaconda中的python，這樣會減小後續代碼實現階段，函數庫的安裝，如numpy,librosa等）
https://blog.csdn.net/u010824101/article/details/85239223 Windows10中用Anaconda的conda環境和VSCode工具來編寫python代碼
https://blog.csdn.net/zzc15806/article/details/79603994 音頻處理庫—librosa的安裝與使用

matlab應用環境：

如何安裝以及破解matlab，請參考下方資料來源（大佬寫的很是詳細），這裏不做贅述了。

安裝前須知（提醒）：

1.安裝全程須斷開網絡，不然破解不成功；

2.解壓和安裝前先關閉360、騰訊管家等殺毒軟件，防止誤殺破解補丁，致使破解失敗；

3.Matlab 2018b適用於WIN7/8/10（64位）系統，親測可用！

4.推薦電腦最低配置：內存8G+，處理器酷睿I5+；

文章參考來源：
https://blog.csdn.net/qq_26900233/article/details/88816789 MATLAB R2018b 安裝教程（含軟件資源）

2、代碼實現

python源碼：

import librosa
import librosa.display
import 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;       % 窗函數長度設置爲128
han_win = hanning(win_sz);        % 選擇海明窗
nfft = 2^nextpow2(length(han_win));       % DFT點數（一般取最接近信號長度的2的整數次冪，即 nfft = 2^nextpow2(length(window)) ）=128=win_size
nooverlap = 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 = 3
c=2*Fc*totalscal;    % 測得c = 1536
scals=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;
figure
imagesc(t,f,abs(coefs));
set(gca,'YDir','normal')
colorbar;
xlabel('時間 t/s');
ylabel('頻率 f/Hz');
title('小波時頻圖')

文章參考來源：
https://blog.csdn.net/lanchunhui/article/details/72240693 matlab 時頻分析（短時傅里葉變換、STFT）
https://blog.csdn.net/weixin_42765703/article/details/104871604 解析MATLAB短時傅里葉變換函數spectrogram()
https://blog.csdn.net/zb1165048017/article/details/80682473 【音頻處理】短時傅里葉變換
https://blog.csdn.net/weixin_44586750/article/details/103219952 數字信號處理——時頻分析（短時傅里葉變換+手動實現時頻分析）
https://blog.csdn.net/qq_42815385/article/details/89095135 STFT（短時傅里葉變換）音頻特徵提取，用於語音識別 python
https://blog.csdn.net/Barry_J/article/details/81065932 音頻特徵提取——python/ librosa工具包使用
https://blog.csdn.net/qq_34218078/article/details/101255636 librosa語音信號處理
https://blog.csdn.net/yunman2012/article/details/103542690 librosa處理音頻信號

3、總結（做者小白話）

寫這篇博客的原因，主要想對近期音頻信號處理相關知識的學習，作一個系統性的梳理，同時也但願可以對你們有所幫助（學習階段，查閱了很多大佬文章，收益匪淺，這些都成爲我寫博客，分享本身的學習歷程，與經驗教訓的一個動力）。
整個學習階段，從瞭解傅里葉變換，到python、matlab環境的配置和使用，以及後續的音頻文件處理，STFT變換等，走過了許多彎路，也遇到了一些不常見的bug（主要是python環境配置，庫安裝時），但最終仍是堅持了下來。目前對音頻處理的相關學習還不是很深刻，文章還存在許多不全面的地方，在後續深刻學習過程當中，會對不合理的地方作出調整，也但願各位大佬能多多指教。相互學習，共同進步啊，哈哈哈。（文章引用了一些大佬的文章（文章小節處有標註），以便各位同仁能更全面的瞭解，固然也會存在一些疏漏，請指教，望見諒。）
補充，寫文章階段時，發現了一位大佬關於python語音基礎操做的文章，內容很詳實，包含音頻生成，語音分析、處理的介紹，收益匪淺，文章最後附原文連接，但願有助於你們理解語音信號處理。

文章參考來源： https://blog.csdn.net/sinat_18131557/article/details/106440757 Python語音基礎操做（對python處理語音信號的過程介紹的很詳細，推薦！）