MATLAB信号处理仿真.docx
《MATLAB信号处理仿真.docx》由会员分享，可在线阅读，更多相关《MATLAB信号处理仿真.docx（5页珍藏版）》请在冰点文库上搜索。
MATLAB信号处理仿真
MATLAB信号处理仿真
MATLAB信号处理仿真
本次，我们将展示另外一个信号分析的有力工具，短时傅里叶变换。
对于探讨经典信号处理的问题，我的感觉相对舒畅，因为从信号与系统到离散时间信号处理，建立在无穷时域前提下的各种结论，完美的仿佛传说中的理想国与乌托邦。
尽管DFT与FFT的时域有限长采样造成了频谱泄露，但是我们似乎仍然心存希望，至少可以通过增加采样长度来降低泄露，至少还有这一途径可以让我们感到些许温暖，对于频率成分固定不变的信号，这种增加采样长度的方法可以让我们逼近理想主义的完美。
但是，对于频率成分随时间发生改变的信号，以上的方法陷入一个两难境地，如果采集很长的一段数据进行DFT分析，可以很精确的得知出现了那些频率成分，但是却无法知道这个成分是在何时出现的，或者，采集很短的数据做DFT分析，可以相对准确的知道在什么时候出现了某个频率成分，然而对于这个成分的具体频率值可能就探测的不那么准确，于是，对于频谱成分随时间变化的信号，我们的傅里叶分析真的变成了一个跷跷板问题，一头好了，另一头就糟糕。
不管怎样，有总比没有强，这就是本次的主题，短时傅里叶变换，在奥本海姆的书中称为“依时
呢？
的确，提高重叠点数是有帮助，但是，得到的分析结果仍然是一种不清晰的结果。
口说无凭，下面我们通过实验现象来观察。
首先我们观察一个DTMF信号离散双音多频信号，说白了就是电话拨号音，这个是我手机号码的电话拨号音，如果你能够译码出我的手机号，欢迎发短信告诉我，请注明是在信号处理博客上看到的，我的博客ID是“duweitao”，请不要把我的手机号公布在网站上，免得给大家添麻烦。
首先看这个拨号音的时域信号，是若干的时域脉冲的组合，脉冲之间是短暂的静默信号。
下图第一行是全局波形，第二行是放大的脉冲局部，每一个拨号数字是用2个频率分量来表示的，不同数字的两个频率分量的组合不同。
然后我们把全部的时域信号拿来做个DFT，我们看到似乎是有6根谱线，但是你能说出那根谱线在何时出现么，或者某个谱线总共出现了几次呢？
接下来，采用256点的短时DFT进行分析，用光谱图spectrogram来观察结果，可以清楚的看到在时间轴上各个断续的音调，但是频率轴上音调的中心位置就不那么清晰了。
然后，我们采用1024点的短时DFT进行分析，这一次可以清楚的看到频率的中心位置，但是，你发现了吗，有些音调之间本来存在的空隙，在这个光谱图上看不见了，可是，事实上它们是存在的。
所以说时域和频域的解析能力是一个跷跷板的两头，你要优化这一头，那么另一头就变得糟糕。
如果你对“spectrogram”函数不了解，请自行在matlab中doc之，或者，请看下面这张图，光谱图是以下这个3D谱图的俯视图，即从上往下看，能量越强，离观察者越近的谱峰，颜色越热。
把同样的方法，应用于附件的“白杨礼赞”普通话播音文件，我们可以看到如下的图片
首先是语谱图
然后是3D视图从此中，我们可以看到，人类的语音中，充满了时变的频谱分量，人类语音信号中有一个重要参数叫做基音周期，亦即基波的频率的倒数，这个基频的频率和其各次倍频的能量关系，以及它们的时变规律，决定了一个人的声音特征。
语言，多么普通，我们天天在使用，用来传递含义或是表达感情。
语音分析与合成是最早使用数字信号处理作为研究工具的领域之一，尽管人类的语音识别、合成的算法日益成熟乃至可以商用，然而，在如此平常的信号里却蕴藏了我们永远也无法完全探究的奥秘，没有两个人的语谱图特征是完全一样的，这东西就像指纹一样，在刑侦学科上，被称作“声纹”，说实话，每次看这种语谱图我的心中都充满敬畏，为什么会是这样？
我不知道，没人知道，也许只有上帝知道。
%///////////////////////////////////////////////////////////
%shorttimeDFTanalysewavfile
%///////////////////////////////////////////////////////////
closeall;
clear;
clc;
%settargetfileatlatterline
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（'Totalsignal'）;
subplot（2,1,2）;
sect_len=800;
plot（wav_data（1:
sect_len））;title（'Firstpart'）;kaiser_win_spectrum_plot（fs,wav_data,kaiser_beta）;
y_min=-160;y_max=10;
ylim（[y_min,y_max]）;%sety-axisrangetitle（'InputWaveFile,NormalizedSpectrum','fontsize',14）;dscp_str=sprintf（'fs%.2fKHz,L%d,q%dbits,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;
%usetictocgetrunningtime
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（'datalen%d\nwinlen%d\noverlap%d\ndftlen%d',...
data_len,len_w,len_o,len_d）;
text（text_x,text_y,dscp_str,...
'FontSize',14,'color','g','BackgroundColor','k'）;
colormap（'jet'）;
%plot3d-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）;
升级会员 