【音频处理】短时傅里叶变换

前言

上一篇博客讲了离散傅里叶变换，里面的实例是对整个信号进行计算，虽然理论上有N点傅里叶变换(本博客就不区分FFT和DFT了，因为它俩就是一个东东，只不过复杂度不同)，但是我个人理解是这个N点是信号前面连续的N个数值，即N点FFT意思就是截取前面N个信号进行FFT，这样就要求我们的前N个采样点必须包含当前信号的一个周期，不然提取的余弦波参数与正确的叠加波的参数相差很大。

如果在N点FFT的时候，如果这N个采样点不包含一个周期呢？或者说我们的信号压根不是一个周期函数咋办？或者有一段是噪音数据呢？如果用FFT计算，就会对整体结果影响很大，然后就有人想通过局部来逼近整体，跟微积分的思想很像，将信号分成一小段一小段，然后对每一小段做FFT，就跟分段函数似的，无数个分段函数能逼近任意的曲线((⊙o⊙)…应该没错吧)，这样每一段都不会互相影响到了。

下面的参考博客中有一篇的一句话很不错：在短时傅里叶变换过程中，窗的长度决定频谱图的时间分辨率和频率分辨率，窗长越长，截取的信号越长，信号越长，傅里叶变换后频率分辨率越高，时间分辨率越差；相反，窗长越短，截取的信号就越短，频率分辨率越差，时间分辨率越好，也就是说短时傅里叶变换中，时间分辨率和频率分辨率之间不能兼得，应该根据具体需求进行取舍。

国际惯例，参考博客：

基于MATLAB短时傅里叶变换和小波变换的时频分析

小波前奏–短时傅里叶变换

短时傅里叶变换原理解

matlab自带的短时傅里叶变换函数spectrogram

理论及实现

其实就是多了几个参数，需要指定的有：

• 每个窗口的长度：nsc
• 每相邻两个窗口的重叠率：nov
• 每个窗口的FFT采样点数：nff

可以计算的有：

• 海明窗：w=hamming(nsc, 'periodic')

• 信号被分成了多少片：$\frac{len(S)-nsc}{nsc-nov}$

• 短时傅里叶变换：
  $$X(k)=\sum _{n=1}^{N}w(n)\ast x(n)\ast exp(-j\ast 2\pi \ast (k-1)\ast (n-1)/N),1<=k<=N$$
  其实和FFT的公式一样，只不过多了个海明窗加权

直接撸代码：

①先设置参数：

%默认设置：
% nsc=floor(L/4.5);%海明窗的长度
% nov=floor(nsc/2);%重叠率
% nff=max(256,2^nextpow2(nsc));%N点采样长度
%也可手动设置
nsc=100;%海明窗的长度,即每个窗口的长度
nov=30;%重叠率
nff=256;%N点采样长度

这里面有个默认设置，就是调用matlab自带的短时傅里叶变换时，如果没指定相关参数，就会采用默认参数值，这个可以去mathwork官网看。

②计算海明窗以及初始化结果值：

h=hamming(nsc, 'periodic');%计算海明窗的数值，给窗口内的信号加权重
coln = 1+fix((L-nsc)/(nsc-nov));%信号被分成了多少个片段
%如果nfft为偶数，则S的行数为(nfft/2+1)，如果nfft为奇数，则行数为(nfft+1)/2
%因为matlab的FFT结果是对称的，只需要一半
rown=nff/2+1;
STFT_X=zeros(rown,coln);%初始化最终结果

这里的信号被划分的片段数目可以按照卷积的方法计算

③对每个片段码公式：

%对每个片段进行fft变换
index=1;%当前片段第一个信号位置在原始信号中的索引
for i=1:coln
    %提取当前片段信号值,并用海明窗进行加权
    temp_S=S(index:index+nsc-1).*h';
    %进行N点FFT变换
    temp_X=fft(temp_S,nff);
    %取一半
    STFT_X(:,i)=temp_X(1:rown)';
    %将索引后移
    index=index+(nsc-nov);
end

可以发现我这里没码公式，因为上一篇博客证明了手撸的DFT与matlab自带的FFT公式一样，有高度强迫症的可以把上一篇博客的DFT写成一个函数，然后把此处的FFT换成你的函数名即可。注意这里的关键操作有两点：

• 对当前窗口的输入信号进行海明加权
• 窗口中输入信号的获取方法有点类似于卷积，卷积核大小是1*nsc，步长是nsc-nov

④正确性验证：与matlab自带的STFT函数spectrogram的结果进行比较：

%% matlab自带函数
[spec_s,spec_f,spec_t]=spectrogram(S,hamming(nsc, 'periodic'),nov,nff,Fs);
%减法，看看差距
plot(abs(spec_s)-abs(STFT_X))

啥也不说了，稳如狗

频谱解读

直接计算每个坐标轴的数值就知道其代表的意思了，其实它和一款叫做Praat的软件所画的图很类似，贴一张图，来源百度

上面是声线，就是直接audioread声音得到的数值，下面就是声谱图，看着是二维图形，其实是3D图形，坐标轴分别代表时间，频率，以及当前时间在当前频率上的幅值。

那么在matlab中如何计算这些值？如下：

回顾一下整个快速离散傅里叶变换的流程：

• 利用窗口和重叠率对整个输入信号进行片段划分
• 对每个片段的信号做N点傅里叶变换，并利用海明窗加权

接下来解析三个坐标，先规定一下横坐标表示频率，纵坐标表示时间，第三个坐标表示幅值：

• 第三个坐标：幅值的计算与FFT的幅值计算一样，都是计算得到的结果除以采样点N，再乘以2，只不过这里要利用海明值进行缩放处理

  % 海明窗口的均值
  K = sum(hamming(nsc, 'periodic'))/nsc;
  %获取幅值(除以采样长度即可，后面再决定那几个除以2)，并依据窗口的海明系数缩放
  STFT_X = abs(STFT_X)/nsc/K;
  % correction of the DC & Nyquist component
  if rem(nff, 2)                     % 如果是奇数，说明第一个是奈奎斯特点，只需对2:end的数据乘以2
      STFT_X(2:end, :) = STFT_X(2:end, :).*2;
  else                                % 如果是偶数，说明首尾是奈奎斯特点，只需对2:end-1的数据乘以2
      STFT_X(2:end-1, :) = STFT_X(2:end-1, :).*2;
  end
  

• 横坐标：频率
  首先要知道当前窗口代表了多少频率，它是在总频率Fs的基础上，对每个窗口进行nff点采样，需要注意的是进行多少次nff采样，当前窗口就有多少个频率值，它们之间是均匀的Fs/nff，这个也通常被称作频率到分辨率的比率。这里看论文《Constant-Q transform toolbox for music processing 》中的一句话：

  The discrete short-time Fourier transform gives a constant resolution for each bin or frequency sampled equal to the sampling rate dividedbythewindowsizein samples. This means,for example,if we takea window of 1024 samples with a sampling rate of 32O30 samples/s (reasonable for musical signals),there solution is 31.3Hz
  

  就是说我们从采样率为32030Hz的样本中挑选包含1024个样本的窗口，那么分辨率就是$\frac{32030}{1024}=31.3Hz$

  所以横坐标的值就是$i\times \frac{Fs}{nff},(i\in (0,nff))$

  %频率轴
  STFT_f=(0:rown-1)*Fs./nff;
  

• 纵坐标时间：
  因为采用了窗口，所以纵坐标的时间比实际时间短，每个坐标代表当前窗口中间信号在原始信号中的索引，窗口是nsc，重叠率是nov，那么每次向前挪的步长为nsc-nov，总共挪coln-1次，需要注意的是这个挪是在采样点上挪，需要将采样点挪转换为时间挪，由于整个信号采样率是Fs，代表每秒的采样数，那么相邻的采样点时间间隔是1/Fs，自然就得到了纵坐标的刻度：

  %时间轴
  STFT_t=(nsc/2:(nsc-nov):nsc/2+(coln-1)*(nsc-nov))/Fs;
  

• 额外信息：赋值转分贝
  我也不清楚这个意义是啥，但是在matlab中有对应函数，而且软件praat和默认函数spectrogram的结果图中都有这个信息，所以还是算一下吧，公式很简单$y=20\ast \log 10(x)$，直接码公式:

  STFT_X = 20*log10(STFT_X + 1e-6);
  

  这里加个很小的值以后，画图好看点

坐标值都得到，直接mesh出来就行，整个代码如下：

%% 画频谱图
% 海明窗口的均值
K = sum(hamming(nsc, 'periodic'))/nsc;
%获取幅值(除以采样长度即可，后面再决定哪几个除以2)，并依据窗口的海明系数缩放
STFT_X = abs(STFT_X)/nsc/K;
% correction of the DC & Nyquist component
if rem(nff, 2)                     % 如果是奇数，说明第一个是奈奎斯特点，只需对2:end的数据乘以2
    STFT_X(2:end, :) = STFT_X(2:end, :).*2;
else                                % 如果是偶数，说明首尾是奈奎斯特点，只需对2:end-1的数据乘以2
    STFT_X(2:end-1, :) = STFT_X(2:end-1, :).*2;
end

% convert amplitude spectrum to dB (min = -120 dB)
%将幅值转换成分贝有函数是ydb=mag2db(y)，这里直接算
STFT_X = 20*log10(STFT_X + 1e-6);
%时间轴
STFT_t=(nsc/2:(nsc-nov):nsc/2+(coln-1)*(nsc-nov))/Fs;
%频率轴
STFT_f=(0:rown-1)*Fs./nff;
% plot the spectrogram
figure
surf(STFT_f, STFT_t,  STFT_X')
shading interp
axis tight
box on
view(0, 90)
set(gca, 'FontName', 'Times New Roman', 'FontSize', 14)
xlabel('Frequency, Hz')
ylabel('Time, s')
% title('Amplitude spectrogram of the signal')
title('手绘图')

handl = colorbar;
set(handl, 'FontName', 'Times New Roman', 'FontSize', 14)
ylabel(handl, 'Magnitude, dB')

对比一下matlab自带函数spectrogram和我们手绘图的正面图和3D图：

从图里面很容易看出来咱们输入的这个音频信号由50和100两个频率组成，幅值大概在10到20，what？咋少了一半，(⊙o⊙)…，后面再研究为啥少了一半还，按理说乘以2了呀，反正频率对了

后记

感觉对于FFT的理解告一段落，先把蝶形算法搁着，下一步就是折腾常Q变换(Constant-Q transform)了，目前的用处是一个音乐的一拍可能有很多音组合而成，但是每个音的频率又不一样，那么就需要设置不同的窗口进行采样，相当于进行了多次STFT操作，只不过每次的窗口大小不同罢了，有兴趣可以看一波论文：《Calculation of a constant Q spectral transform 》，有张图介绍了CQT和DFT的区别，具体我还在研究，感觉就是把STFT的for循环里面的N变成动态计算的就行了。

贴一波代码，直接贴了，懒得贴网盘：

%短时傅里叶变换STFT
%依据FFT手动实现STFT
clear
clc
close all
Fs = 1000;            % Sampling frequency
T = 1/Fs;             % Sampling period
L = 1000;             % Length of signal
t = (0:L-1)*T;        % Time vector
S = 20*cos(100*2*pi*t)+40*cos(50*2*pi*t);%0.2-0.7*cos(2*pi*50*t+20/180*pi) + 0.2*cos(2*pi*100*t+70/180*pi) ;

%% 所需参数
%主要包含:信号分割长度(默认分割8个窗口)，海明窗口，重叠率，N点采样
%默认设置：
% nsc=floor(L/4.5);%海明窗的长度
% nov=floor(nsc/2);%重叠率
% nff=max(256,2^nextpow2(nsc));%N点采样长度
%也可手动设置
nsc=100;%海明窗的长度,即每个窗口的长度
nov=0;%重叠率
nff=max(256,2^nextpow2(nsc));%N点采样长度
%% 手动实现STFT
h=hamming(nsc, 'periodic');%计算海明窗的数值，给窗口内的信号加权重
coln = 1+fix((L-nsc)/(nsc-nov));%信号被分成了多少个片段
%如果nfft为偶数，则S的行数为(nfft/2+1)，如果nfft为奇数，则行数为(nfft+1)/2
%因为matlab的FFT结果是对称的，只需要一半
rown=nff/2+1;
STFT_X=zeros(rown,coln);%初始化最终结果
%对每个片段进行fft变换
index=1;%当前片段第一个信号位置在原始信号中的索引
for i=1:coln
    %提取当前片段信号值,并用海明窗进行加权
    temp_S=S(index:index+nsc-1).*h';
    %进行N点FFT变换
    temp_X=fft(temp_S,nff);
    %取一半
    STFT_X(:,i)=temp_X(1:rown)';
    %将索引后移
    index=index+(nsc-nov);
end

%% matlab自带函数
spectrogram(S,hamming(nsc, 'periodic'),nov,nff,Fs);
title('spectrogram函数画图')
[spec_X,spec_f,spec_t]=spectrogram(S,hamming(nsc, 'periodic'),nov,nff,Fs);
%减法，看看差距
% plot(abs(spec_X)-abs(STFT_X))

%% 画频谱图
% 海明窗口的均值
K = sum(hamming(nsc, 'periodic'))/nsc;
%获取幅值(除以采样长度即可，后面再决定哪几个乘以2)，并依据窗口的海明系数缩放
STFT_X = abs(STFT_X)/nsc/K;
% correction of the DC & Nyquist component
if rem(nff, 2)                     % 如果是奇数，说明第一个是奈奎斯特点，只需对2:end的数据乘以2
    STFT_X(2:end, :) = STFT_X(2:end, :).*2;
else                                % 如果是偶数，说明首尾是奈奎斯特点，只需对2:end-1的数据乘以2
    STFT_X(2:end-1, :) = STFT_X(2:end-1, :).*2;
end

% convert amplitude spectrum to dB (min = -120 dB)
%将幅值转换成分贝有函数是ydb=mag2db(y)，这里直接算
STFT_X = 20*log10(STFT_X + 1e-6);
%时间轴
STFT_t=(nsc/2:(nsc-nov):nsc/2+(coln-1)*(nsc-nov))/Fs;
%频率轴
STFT_f=(0:rown-1)*Fs./nff;
% plot the spectrogram
figure
surf(STFT_f, STFT_t,  STFT_X')
shading interp
axis tight
box on
view(0, 90)
set(gca, 'FontName', 'Times New Roman', 'FontSize', 14)
xlabel('Frequency, Hz')
ylabel('Time, s')
title('Amplitude spectrogram of the signal')

handl = colorbar;
set(handl, 'FontName', 'Times New Roman', 'FontSize', 14)
ylabel(handl, 'Magnitude, dB')

博客已同步更新到微信公众号中，有问题可以在微信公众号中私聊，或者在csdn博客下面评论。