必威体育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)nscnscnov\frac{len(S)-nsc}{nsc-nov}nsc−novlen(S)−nsc​

  • 短时傅里叶变换:

    X(k)=n=1Nw(n)x(n)exp(j2π(k1)(n1)/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=20log10(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值,虽然百度权重只是爱站网提供的一个评

分享到:

栏目导航

推荐阅读

热门阅读