MATLAB信号处理仿真.doc

MATLAB信号处理仿真 本次，我们将展示另外一个信号分析的有力工具，短时傅里叶变换。 对于探讨经典信号处理的问题，我的感觉相对舒畅，因为从信号与系统到离散时间信号处理，建立在无穷时域前提下的各种结论，完美的仿佛传说中的理想国与乌托邦。尽管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 linewav_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 rangetitle(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 timetic ;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 %dn overlap %dn dft len %d,. data_len, len_w, len_o, len_d);text(text_x, text_y,dscp_str, . FontSize,14, color, g,BackgroundColor,k) ;co