必威体育Betway必威体育官网
当前位置：首页 > IT技术

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

时间:2019-06-25 16:44:14来源:IT技术作者:seo实验室小编阅读:66次「手机版」
　

短时傅里叶变换

前言

上一篇博客讲了离散傅里叶变换，里面的实例是对整个信号进行计算，虽然理论上有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')

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

  • 短时傅里叶变换：

    X(k)=∑n=1Nw(n)∗x(n)∗exp(−j∗2π∗(k−1)∗(n−1)/N),1&lt;=k&lt;=N X(k) = \sum_{n=1}^N w(n)* x(n)*exp(-j*2\pi*(k-1)*(n-1)/N), 1 &lt;= k &lt;= N X(k)=n=1∑N​w(n)∗x(n)∗exp(−j∗2π∗(k−1)∗(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 pidedbythewindowsizein 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个样本的窗口，那么分辨率就是320301024=31.3Hz\frac{32030}{1024}=31.3Hz102432030​=31.3Hz

    所以横坐标的值就是i×Fsnff,(i∈(0,nff))i\times \frac{Fs}{nff},(i\in(0,nff))i×nffFs​,(i∈(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∗log⁡10(x)y=20*\log10(x)y=20∗log10(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')

相关阅读

时频分析：短时傅里叶变换实现（2）

目录： 前言 实验环境 Matlab spectrogram函数 1语法 2使用说明 3代码如下： 3.1坐标轴变换 3.2鸟声信号的功率谱密度 3.2.1平方鸟声

快速傅里叶变换及java实现

傅里叶变换：傅里叶变换是一种线性的积分变换。它的理论依据是：任何连续周期信号都可以由一组适当的正弦曲线组合而成，即使用简单的正

傅里叶变换及其实现（MATLAB）

傅立叶变换 傅立叶变换是一种常见的分析方法，傅立叶变换将满足一定条件的函数表示为一些函数的加权和（或者积分）。可以分为四个类别

实战教你如何在短时间打造一个男性流量站

说到流量站，有很多种，有做电影的，有做图片的，还有做种子站啥的，很多很多，就不一一罗列了，一般电影站和图片站相对SEO来说会难一点，因为大

SEO技巧：如何在短时间内提升百度权重？

外国人用谷歌、中国用百度，做百度SEO的站长绝大部分非常在意网站的权重 ，也就是我们所说的BR值，虽然百度权重只是爱站网提供的一个评

赞()
分享到:

猜你喜欢

栏目导航

推荐阅读

热门阅读

seo实验室提供seo优化服务，有网站优化、网站建设、运营推广等需求请联系我们!
Copyright © 2001-2017 .seo实验室 版权所有 网站地图