Matlab_spectrogram_短时傅里叶分析_实现与讨论 [未完成]
Posted adgk07
tags:
篇首语:本文由小常识网(cha138.com)小编为大家整理,主要介绍了Matlab_spectrogram_短时傅里叶分析_实现与讨论 [未完成]相关的知识,希望对你有一定的参考价值。
在语音与音乐处理过程中,常用到短时傅里叶变换(Short Time Fourier Transformation, STFT)。在一些学习路径中,STFT也是学习小波之前的预备知识。本文简单实现了 Matlab 中 Spectrogram 函数,在没有小波知识支撑下,讨论了参数的选择,以及分辨率相关的问题。
参考博客:
来源: 作者:
来源: 作者:
来源: 作者:
短时傅里叶变换简介
自己写一下对STFT的认识,不一定对。
短时傅里叶变换,其实还是傅里叶变换,只不过把一段长信号按信号长度(nsc)、重叠点数(nov)重新采样。原始信号的每一个数据点,都有一个电压值对应,即只有一个直流自由度。重新采样之后,数据点个数变少,每一个数据点由(nsc)个点组成,即得到了其他频率的自由度。但是这个点的频谱是有上下限的——上限受采样频率界定,下限受数据时长界定,和 FFT 一样。
这么处理数据是有缘由的。语音与音乐信号中,信号频率常常变化,而且频率成分丰富,导致简单的傅里叶变换不能很好的描述信号。然而我们人耳处理声音的能力很强,可以对很短的一段声音进行精确的频率分析。所以在语音识别以及语音信号预处理过程中,STFT 是 FFT 的仿生改良版(个人理解)。当然,在其他声学相关方向中,STFT也是蛮有用的。没查过相关文献不敢乱说,模态分析、瞬态特性等应该能用上(个人猜测)。
上图是我用本文的例程实现的语音“XX大学”的 STFT 结果。语音是自己拿录音笔录的,时间域上没有做截断,但由于自己声线太低,频率域截到 2000 Hz。可以看到,我在开始之后等了一段时间才开始说话,后面有几段“嘎嘣脆”的噪声。还可以看到,我的元音发音还是很清晰的,但声线真心低。理论上能从这张图还原出原始信号,不在话下。
Spectrogram 函数函数实现
简单编了一个 STFT 函数如下:
function [ S , W , T ] = mf_spectrogram... ( signal , nsc , nov , fs ) %MF_SPECTROGRAM Short-time Fourier transform realization % Input % signal - Signal vector % nsc - Abb. of Number SeCtion % nov - Abb. of Number OverLap % fs - Abb. of Frequency of Sample % Output % S - A matrix that each colum is a FFT for time of nsc % W - A vector labeling frequency % T - A vector labeling time % Signal Preprocessing h = hamming(nsc, ‘periodic‘); % Hamming weight function L = length(signal); % Length of Signal nst = nsc-nov; % Number of STep per colum ncl = fix( (L-nsc)/nst ) + 1; % Number of CoLum nff = 2^nextpow2(nsc); % Number of Fast Fourier Transformation Ps = zeros(nsc,ncl); for n = 1:ncl % Ps means Processed Signal Ps(:,n) = signal( (n-1)*nst+1 : (n-1)*nst+nsc ).*h‘; end % Ps is a matrix % Short-time Fourier transform STFT0 = fft(Ps,nff); % Turn into DFT in dB STFT1 = abs(STFT0/nff); STFT2 = STFT1(1:nff/2+1,:); % Symmetric results of FFT STFT2(2:end-1,:) = 2*STFT2(2:end-1,:); % Add the value of the other half STFT3 = 20*log10(STFT2); % Turn sound pressure into dB level % Axis Generating fax = fs*(0:(nff/2))/nff; % Frequency axis setting tax = ( .5*nsc : nst : nst*(ncl-1)+.5*nsc ) / fs; % Time axis generating [ffax,ttax] = meshgrid(tax,fax); % Generating grid of figure % Output W = ffax; T = ttax; S = STFT3; end
为了节省代码量,我搞了一个绘图函数:
function [ ] = my_pcolor( f , t , s ) %MY_PCOLOR 绘图函数 % Input % f - 频率轴矩阵 % t - 时间轴矩阵 % s - 幅度轴矩阵 gca = pcolor(f,t,s); % 绘制伪彩色图 set(gca, ‘LineStyle‘,‘none‘); % 取消网格,避免一片黑 handl = colorbar; % 彩图坐标尺 set(handl, ‘FontName‘, ‘Times New Roman‘, ‘FontSize‘, 14) ylabel(handl, ‘Magnitude, dB‘) % 坐标尺名称 end
下面是实现的例程:
clc clear close all % 基本参数 fa = [ 50 800 ]; % 扫描频率上下限 fs = 6400; % 采样频率 % 分辨率相关设定参数 te = 1; % [s] 扫频时间长度 fre = 8; % [s] 频率分辨率 tre = 0.002; % [Hz] 时间分辨率 % Chirp 信号生成 t = 0:1/fs:te; % [s] 时间序列 sc = chirp(t,fa(1),te,fa(2)); % 信号生成 % 分辨率相关输入参数 nsc = floor(fs/fre); % nov = floor(nsc-(fs*tre)); nov = floor(nsc*0.9); nff = max(256,2^nextpow2(nsc)); % 计算与绘制结果 subplot(1,3,1) % 绘制自编函数结果 [S,W,T] = mf_spectrogram(sc,nsc,nov,fs); my_pcolor(W,T,S) caxis([-130.86,-13.667]); title(‘自编函数‘); xlabel(‘时间 second‘); ylabel(‘频率 Hz‘); subplot(1,3,2) % 绘制 Matlab 函数结果 s = spectrogram(sc,hamming(nsc),nov,nff,fs,‘yaxis‘); % Turn into DFT in dB s1 = abs(s/nff); s2 = s1(1:nff/2+1,:); % Symmetric results of FFT s2(2:end-1,:) = 2*s2(2:end-1,:); % Add the value of the other half s3 = 20*log10(s2); % Turn sound pressure into dB level my_pcolor(W,T,s3 ) caxis([-130.86,-13.667]); title(‘Matlab 自带函数‘); xlabel(‘时间 second‘); ylabel(‘频率 Hz‘); subplot(1,3,3) % 两者误差 my_pcolor(W,T,20*log10(abs(10.^(s3/20)-10.^(S/20)))) caxis([-180,-13.667]); title(‘error‘); ylabel(‘频率 Hz‘); xlabel(‘时间 second‘); suptitle(‘Spectrogram 实现与比较‘);
跑的结果我就不贴了,Demo确定是没问题的。网上也有很多相关话题的,例程都比较简单,但我非常善于把问题复杂化:te(扫描长度)、fre(频率下限)、tre(时域分辨率)、nsc(单段数据长度)、nov(重叠数据点数)、nff(FFT点数) 等参数都是为了看设定什么样的参数能够使得 STFT 结果最优。我只是科普性地了解过小波,公式推导还没展开,所以在此只能“实验性”地讨论“尺度”相关的话题。
观察频率上限受到采样频率限制,频率下限受到nsc限制。nsc越大,单段数据时间跨度越长,一旦频率变化率较快,结果就很粗(te = 1 ; fre = 2 ; nov = 99%)。另一方面,为了保持 FFT 频率分辨率,设置 nff 不低于256,这使得当nsc较小时,单段信号中,低频成分可能也就几个周期就结束了,造成了频域上的失真。一句话,单段信号过短,低频效果不好;单段信号过长,捕捉不到频率快速变化的信号。
FFT 周期数、点数、补零数对结果的影响
以上是关于Matlab_spectrogram_短时傅里叶分析_实现与讨论 [未完成]的主要内容,如果未能解决你的问题,请参考以下文章
短时傅里叶变换(Short Time Fourier Transform)原理及 Python 实现