为什么要做手机网站,wordpress 获取插件目录,高端企业网站开发,做品牌推广应该怎么做前言
上一篇博客讲了离散傅里叶变换#xff0c;里面的实例是对整个信号进行计算#xff0c;虽然理论上有N点傅里叶变换(本博客就不区分FFT和DFT了#xff0c;因为它俩就是一个东东#xff0c;只不过复杂度不同)#xff0c;但是我个人理解是这个N点是信号前面连续的N个数值…前言
上一篇博客讲了离散傅里叶变换里面的实例是对整个信号进行计算虽然理论上有N点傅里叶变换(本博客就不区分FFT和DFT了因为它俩就是一个东东只不过复杂度不同)但是我个人理解是这个N点是信号前面连续的N个数值即N点FFT意思就是截取前面N个信号进行FFT这样就要求我们的前N个采样点必须包含当前信号的一个周期不然提取的余弦波参数与正确的叠加波的参数相差很大。
如果在N点FFT的时候如果这N个采样点不包含一个周期呢或者说我们的信号压根不是一个周期函数咋办或者有一段是噪音数据呢如果用FFT计算就会对整体结果影响很大然后就有人想通过局部来逼近整体跟微积分的思想很像将信号分成一小段一小段然后对每一小段做FFT就跟分段函数似的无数个分段函数能逼近任意的曲线((⊙o⊙)…应该没错吧)这样每一段都不会互相影响到了。
下面的参考博客中有一篇的一句话很不错在短时傅里叶变换过程中窗的长度决定频谱图的时间分辨率和频率分辨率窗长越长截取的信号越长信号越长傅里叶变换后频率分辨率越高时间分辨率越差相反窗长越短截取的信号就越短频率分辨率越差时间分辨率越好也就是说短时傅里叶变换中时间分辨率和频率分辨率之间不能兼得应该根据具体需求进行取舍。
国际惯例参考博客
基于MATLAB短时傅里叶变换和小波变换的时频分析
小波前奏–短时傅里叶变换
短时傅里叶变换原理解
matlab自带的短时傅里叶变换函数spectrogram
理论及实现
其实就是多了几个参数需要指定的有
每个窗口的长度nsc每相邻两个窗口的重叠率nov每个窗口的FFT采样点数nff
可以计算的有 海明窗whamming(nsc, periodic) 信号被分成了多少片len(S)−nscnsc−nov\frac{len(S)-nsc}{nsc-nov}nsc−novlen(S)−nsc 短时傅里叶变换 X(k)∑n1Nw(n)∗x(n)∗exp(−j∗2π∗(k−1)∗(n−1)/N),1kNX(k) \sum_{n1}^N w(n)* x(n)*exp(-j*2\pi*(k-1)*(n-1)/N), 1 k N X(k)n1∑Nw(n)∗x(n)∗exp(−j∗2π∗(k−1)∗(n−1)/N),1kN 其实和FFT的公式一样只不过多了个海明窗加权
直接撸代码
①先设置参数
%默认设置
% nscfloor(L/4.5);%海明窗的长度
% novfloor(nsc/2);%重叠率
% nffmax(256,2^nextpow2(nsc));%N点采样长度
%也可手动设置
nsc100;%海明窗的长度,即每个窗口的长度
nov30;%重叠率
nff256;%N点采样长度这里面有个默认设置就是调用matlab自带的短时傅里叶变换时如果没指定相关参数就会采用默认参数值这个可以去mathwork官网看。
②计算海明窗以及初始化结果值
hhamming(nsc, periodic);%计算海明窗的数值给窗口内的信号加权重
coln 1fix((L-nsc)/(nsc-nov));%信号被分成了多少个片段
%如果nfft为偶数则S的行数为(nfft/21)如果nfft为奇数则行数为(nfft1)/2
%因为matlab的FFT结果是对称的只需要一半
rownnff/21;
STFT_Xzeros(rown,coln);%初始化最终结果这里的信号被划分的片段数目可以按照卷积的方法计算
③对每个片段码公式
%对每个片段进行fft变换
index1;%当前片段第一个信号位置在原始信号中的索引
for i1:coln%提取当前片段信号值,并用海明窗进行加权temp_SS(index:indexnsc-1).*h;%进行N点FFT变换temp_Xfft(temp_S,nff);%取一半STFT_X(:,i)temp_X(1:rown);%将索引后移indexindex(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的数据乘以2STFT_X(2:end, :) STFT_X(2:end, :).*2;
else % 如果是偶数说明首尾是奈奎斯特点只需对2:end-1的数据乘以2STFT_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个样本的窗口那么分辨率就是32030102431.3Hz\frac{32030}{1024}31.3Hz10243203031.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的结果图中都有这个信息所以还是算一下吧公式很简单y20∗log10(x)y20*\log10(x)y20∗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的数据乘以2STFT_X(2:end, :) STFT_X(2:end, :).*2;
else % 如果是偶数说明首尾是奈奎斯特点只需对2:end-1的数据乘以2STFT_X(2:end-1, :) STFT_X(2:end-1, :).*2;
end% convert amplitude spectrum to dB (min -120 dB)
%将幅值转换成分贝有函数是ydbmag2db(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到20what咋少了一半(⊙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*t20/180*pi) 0.2*cos(2*pi*100*t70/180*pi) ;%% 所需参数
%主要包含:信号分割长度(默认分割8个窗口)海明窗口重叠率N点采样
%默认设置
% nscfloor(L/4.5);%海明窗的长度
% novfloor(nsc/2);%重叠率
% nffmax(256,2^nextpow2(nsc));%N点采样长度
%也可手动设置
nsc100;%海明窗的长度,即每个窗口的长度
nov0;%重叠率
nffmax(256,2^nextpow2(nsc));%N点采样长度
%% 手动实现STFT
hhamming(nsc, periodic);%计算海明窗的数值给窗口内的信号加权重
coln 1fix((L-nsc)/(nsc-nov));%信号被分成了多少个片段
%如果nfft为偶数则S的行数为(nfft/21)如果nfft为奇数则行数为(nfft1)/2
%因为matlab的FFT结果是对称的只需要一半
rownnff/21;
STFT_Xzeros(rown,coln);%初始化最终结果
%对每个片段进行fft变换
index1;%当前片段第一个信号位置在原始信号中的索引
for i1:coln%提取当前片段信号值,并用海明窗进行加权temp_SS(index:indexnsc-1).*h;%进行N点FFT变换temp_Xfft(temp_S,nff);%取一半STFT_X(:,i)temp_X(1:rown);%将索引后移indexindex(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的数据乘以2STFT_X(2:end, :) STFT_X(2:end, :).*2;
else % 如果是偶数说明首尾是奈奎斯特点只需对2:end-1的数据乘以2STFT_X(2:end-1, :) STFT_X(2:end-1, :).*2;
end% convert amplitude spectrum to dB (min -120 dB)
%将幅值转换成分贝有函数是ydbmag2db(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博客下面评论。
