MATLAB信号处理仿真-短时傅里叶变换,时、频分辨力的跷跷板

EDN China2018-10-29 19:37:19

关注我们随时获精品加关注

本次,我们将展示另外一个信号分析的有力工具,短时傅里叶变换


对于探讨经典信号处理的问题,我的感觉相对舒畅,因为从信号与系统到离散时间信号处理,建立在无穷时域前提下的各种结论,完美的仿佛传说中的理想国与乌托 邦。尽管DFT与FFT的时域有限长采样造成了频谱泄露,但是我们似乎仍然心存希望,至少可以通过增加采样长度来降低泄露,至少还有这一途径可以让我们感 到些许温暖,对于频率成分固定不变的信号,这种增加采样长度的方法可以让我们逼近理想主义的完美。但是,对于频率成分随时间发生改变的信号,以上的方法陷 入一个两难境地,如果采集很长的一段数据进行DFT分析,可以很精确的得知出现了那些频率成分,但是却无法知道这个成分是在何时出现的,或者,采集很短的 数据做DFT分析,可以相对准确的知道在什么时候出现了某个频率成分,然而对于这个成分的具体频率值可能就探测的不那么准确,于是,对于频谱成分随时间变 化的信号,我们的傅里叶分析真的变成了一个跷跷板问题,一头好了,另一头就糟糕。不管怎样,有总比没有强,这就是本次的主题,短时傅里叶变换,在奥本海姆 的书中称为“依时傅里叶变换”。 具体的数学推导此处就不再赘述了,并且没有数学排版工具也不方便展示公式,维基百科英文版的短时傅里叶变换写的也比较生动,此处我们定性的讲一讲,就是如 果有一个很长的时间信号序列,我们每次从中取出一个长度为L的片段,加窗做DFT,看看频谱,这就是短时傅里叶变换,由于多引入了一个时间的维度,所以存 在一个问题,就是两次做DFT的L点数据帧,它们在原始信号中的相对位置是怎样的,是否允许有些重叠呢?我们思考一下,如果两帧信号完全不重叠,那么由于 加窗的原因(想想,窗函数的边界处通常都是接近0的数据)交界处的信号突变会被窗函数抹杀掉,所以通常来说两帧数据是要有重叠的,而且重叠的部分越多,对 于捕获突发的信号越有帮助,这样做的代价呢,就是计算量的增加。

那么可能您会有疑问,在做短时傅里叶变换的时候,增加DFT的帧长度L会不会是个一本万利的事情呢?很遗憾,不是,因为对于L点的DFT的分析中找到的谱 峰,其有效意义对应于整个的L点数据,所以无法分辨出是在L点数据中的那个位置出现的,也许你会想,可以把两帧数据的重叠点数增加啊,这样是否利用观测某 个谱峰在那个时间出现呢?的确,提高重叠点数是有帮助,但是,得到的分析结果仍然是一种不清晰的结果。


口说无凭,下面我们通过实验现象来观察。 首先我们观察一个DTMF信号离散双音多频信号,说白了就是电话拨号音,这个是我手机号码的电话拨号音,如果你能够译码出我的手机号,欢迎发短信告诉我,请注明是在信号处理博客上看到的,我的博客ID是“duweitao”,请不要把我的手机号公布在网站上,免得给大家添麻烦。


首先看这个拨号音的时域信号,是若干的时域脉冲的组合,脉冲之间是短暂的静默信号。下图第一行是全局波形,第二行是放大的脉冲局部,每一个拨号数字是用2个频率分量来表示的,不同数字的两个频率分量的组合不同。

、然后我们把全部的时域信号拿来做个DFT,我们看到似乎是有6根谱线,但是你能说出那根谱线在何时出现么,或者某个谱线总共出现了几次呢?


接下来,采用256点的短时DFT进行分析,用光谱图spectrogram来观察结果,可以清楚的看到在时间轴上各个断续的音调,但是频率轴上音调的中心位置就不那么清晰了。


然后,我们采用1024点的短时DFT进行分析,这一次可以清楚的看到频率的中心位置,但是,你发现了吗,有些音调之间本来存在的空隙,在这个光谱图上看不见了,可是,事实上它们是存在的。


所以说时域和频域的解析能力是一个跷跷板的两头,你要优化这一头,那么另一头就变得糟糕。如果你对“spectrogram”函数不了解,请自行在 matlab中doc 之,或者,请看下面这张图,光谱图是以下这个3D谱图的俯视图,即从上往下看,能量越强,离观察者越近的谱峰,颜色越热。


把同样的方法,应用于附件的“白杨礼赞”普通话播音文件,我们可以看到如下的图片首先是语谱图。


然后是3D视图


从此中,我们可以看到,人类的语音中,充满了时变的频谱分量,人类语音信号中有一个重要参数叫做基音周期,亦即基波的频率的倒数,这个基频的频率和其各次倍频的能量关系,以及它们的时变规律,决定了一个人的声音特征。 语言,多么普通,我们天天在使用,用来传递含义或是表达感情。语音分析与合成是最早使用数字信号处理作为研究工具的领域之一,尽管人类的语音识别、合成的算法日益成熟乃至可以商用,然而,在如此平常的信号里却蕴藏了我们永远也无法完全探究的奥秘,没有两个人的语谱图特征是完全一样的,这东西就像指纹一样,在刑侦学科上,被称作“声纹”,说实话,每次看这种语谱图我的心中都充满敬畏,为什么会是这样? 我不知道,没人知道,也许只有上帝知道。

%///////////////////////////////////////////////////////////
% short time DFT analyse wav file
%///////////////////////////////////////////////////////////
close all;
clear;
clc;
% set target file at latter line
wav_fname = 'baiyang_test.wav' ;
wav_fname = 'dial_tone.wav' ;

[wav_data, fs, q_bits] = wavread(wav_fname) ;
data_len = length(wav_data) ;
kaiser_beta = 8 ;

stdft_len_win = 256 ;
stdft_len_overlap = stdft_len_win-16 ;
stdft_len_dft = stdft_len_win ;

figure;
subplot(2,1,1);
plot(wav_data); title('Total signal');
subplot(2,1,2);
sect_len = 800;
plot(wav_data(1:sect_len));title('First part');

kaiser_win_spectrum_plot(fs, wav_data, kaiser_beta);
y_min = -160; y_max = 10;
ylim([y_min, y_max]) ; % set y-axis range

title('Input Wave File, Normalized Spectrum', 'fontsize', 14);

dscp_str = sprintf('fs %.2f KHz, L %d, q %d bits, beta %d',...
fs/1E3, data_len, q_bits, kaiser_beta);

text(0.5*(-fs/2),0.8*y_min,dscp_str, ...
'FontSize',14, 'color', 'g','BackgroundColor','k') ;

len_w = stdft_len_win ;
len_o = stdft_len_overlap ;
len_d = stdft_len_dft ;
figure ;
% use tic toc get running time
tic ;
spectrogram(wav_data, len_w, len_o, len_d, fs);
title('Spectrogram', 'fontsize', 14);
toc ;
text_x = 0.7*(fs/2);
text_y = 0.5* (data_len/fs);
dscp_str = sprintf(' data len %d \n win len %d\n overlap %d\n dft len %d',...
data_len, len_w, len_o, len_d);
text(text_x, text_y,dscp_str, ...
'FontSize',14, 'color', 'g','BackgroundColor','k') ;
colormap('jet');
% plot 3d-mesh
S = spectrogram(wav_data, len_w, len_o, len_d, fs);
figure; mesh(log10(1E-4+abs(S))*10);
xlabel('Time', 'FontSize',14);
ylabel('Frequency', 'FontSize',14);


附件: Matlab_Code_STDFT_WAV.rar请点击阅读原文下载。


回复:傅里叶、阻抗、面试、电源、FPGA、CAN 查看更多好文

喜欢
分享收藏
or



点击下方“阅读原文”更有用
Copyright © 温县电话机虚拟社区@2017